Eric Ferrage

Bruno Lanson
email: bruno.lanson@obs.ujf-grenoble.fr

Natalie Malikova

Alain Plançon

Boris A Sakharov

Victor A Drits

New insights on the distribution of interlayer water in bi-hydrated smectite from X-ray diffraction profile modeling of 00ℓ reflections

come

New insights on the distribution of interlayer water in

bi-hydrated smectite from X-ray diffraction profile modeling 00l reflections

Introduction

Smectite is a 2:1 phyllosilicate whose layer structure consists of an octahedral sheet sandwiched in-between two siliceous tetrahedral sheets. Isomorphic substitutions in either tetrahedral or octahedral sites induce a permanent negative layer charge, which is compensated for by the presence of hydrated cations in the interlayer. The observation of 00ℓ basal reflections on X-ray diffraction (XRD) patterns has shown that with increasing relative humidity smectite expands stepwise, the different steps corresponding to the intercalation of 0, 1, 2 or 3 sheets of H 2 O molecules in the interlayer. [START_REF] Nagelschmidt | [END_REF][START_REF] Bradley | [END_REF][START_REF] Mooney | [END_REF][5][6] From these pioneer studies, it is now commonly accepted that the expandability of 2:1 phyllosilicates is controlled by factors such as the nature of interlayer cations, and the layer charge and its location (octahedral vs. tetrahedral). These general observations have led to different models in which crystalline swelling is controlled by the balance between the repulsive forces between neighboring 2:1 layers and the attractive forces between hydrated interlayer cations and the negatively-charged surface of siloxane sheets. [6][7][8][9][10][11] The development of XRD modeling techniques allowed investigating structures in which different hydration states coexist thus improving these early observations. [START_REF] Ben Brahim | In 5th Meeting of the European Clay Groups[END_REF][START_REF] Ben Brahim | [END_REF][14][15][16][17] Ferrage et al. used such a modeling approach to characterize the hydration of several montmorillonite and beidellite samples and observed that the nature of the interlayer cation, and in particular its affinity for water, influences the layer thickness of bi-hydrated and monohydrated layers. 18,19 They also confirmed that the relative proportions of the different layer types, which correspond to the different hydration states, depend on both the amount and the location of smectite layer charge. In addition, these authors showed that XRD peak profiles and position can be satisfactorily reproduced, especially over the low-angle region (~5-12°2θ Cu Kα), only if hydration heterogeneity is taken into account. They were thus able to refine the structure of smectite and in particular to investigate atomic positions of interlayer species.

In particular, they showed that the atomic positions reported by Moore and Reynolds for H 2 O molecules in bi-hydrated layers induce a dramatic misfit over the medium-to high-angle region (12-50°2θ Cu Kα) by strongly modifying the intensity ratio between the different 00ℓ reflections. [START_REF] Moore | Jr X-ray Diffraction and the Identification and Analysis of Clay Minerals[END_REF]18 The present article thus aims at refining further the structure of interlayer H 2 O in bihydrated smectites from the fit of experimental XRD patterns. The proposed structure is compared with the positional distribution commonly derived from Monte-Carlo simulations, whereas the adjusted amounts of interlayer water are compared with those determined experimentally from water vapor adsorption-desorption experiments.

Background

Smectite hydration heterogeneity as seen by XRD profile modeling. In agreement with the stepwise evolution of the d 001 basal spacing on XRD patterns, the hydration state of smectite has been described using three layer types exhibiting different layer thickness corresponding to the common hydration states reported for montmorillonite in non-saturated conditions. Dehydrated layers (0W -Layer thickness ~9.6-10.1 Å), mono-hydrated layers (1W -Layer thickness ~12.3-12.7 Å), and bi-hydrated layers (2W -Layer thickness ~15.1-15.8 Å) have thus been defined. In the first two layer types, interlayer cations are located in the mid-plane of the interlayer, together with H 2 O molecules for 1W layers. For 2W layers, interlayer cations are also commonly assumed to be located in the mid-plane of the interlayer. [START_REF] Moore | Jr X-ray Diffraction and the Identification and Analysis of Clay Minerals[END_REF] In addition, it is usually assumed that two planes of H 2 O molecules, each bearing 0.69 H 2 O per O 20 (OH) 4 , are located at 0.35 and 1.06 Å from the cation along the c* axis (Debye-Waller parameter B wat ~2 Å 2 for these two planes), whereas a third denser plane (1.20 H 2 O per O 20 (OH) 4 ) is located further from the central interlayer cation at 1.20 Å along the c* axis (B wat = 11 Å 2 ). [START_REF] Moore | Jr X-ray Diffraction and the Identification and Analysis of Clay Minerals[END_REF] The pattern calculated for the Ca-saturated reference SWy-1 montmorillonite (Ca-SWy-1) assuming a homogeneous 2W hydration state and the above configuration for interlayer species is compared on Figure 1a to the experimental pattern recorded at 80% RH. With these usual hypotheses, the calculated pattern fits most of the experimental pattern features but significant discrepancies can be observed over the mediumto the high-angle region in spite of the low intensity diffracted. In particular, the position of the 005 reflection and the low-angle "tail" of the 002 reflection are not well reproduced (Figure 1a). Ferrage et al. challenged this usual configuration of interlayer species, and proposed an alternative configuration that includes a unique plane of H 2 O molecules located at 1.20 Å, along the c* axis, on either side of the central interlayer cation (2WS configuration). 18 The use of this 2WS configuration helps reducing the discrepancies observed for the 003-005 reflections. In particular, this configuration allows decreasing the relative intensity of the 003 and 004 reflections whereas the intensity of the 005 one is increased (Figure 1b). However, in the high-angle region the intensity ratio between the 007 and 008 reflections measured on the calculated pattern is inconsistent with that determined experimentally, although the intensity of the 008 reflection is correctly reproduced.

Ferrage et al. also demonstrated that the common hypothesis of a homogeneous hydration state for smectite is not consistent with the likely existence in smectite of structural heterogeneities affecting the layer charge distribution (from one interlayer to the other or within a given interlayer) and/or location (octahedral vs. tetrahedral). 18 In turn these heterogeneities lead to the coexistence of different layer types in a single structure. Such hydration heterogeneity has been evidenced from the profile modeling of XRD patterns recorded on hydrated smectites. [START_REF] Ben Brahim | In 5th Meeting of the European Clay Groups[END_REF][START_REF] Ben Brahim | [END_REF][14][15][16][17] Ferrage et al. have shown that this heterogeneity is systematically observed whatever the interlayer cation, the relative humidity (RH), and the amount and location of the layer charge deficit. 18,19 It is thus essential to account for the hydration heterogeneity to satisfactorily reproduce the experimental positions and profiles of reflections.

Ferrage et al. have shown indeed that accounting for smectite hydration heterogeneity allows fitting better the profiles of all experimental 00ℓ reflections. 18,19 In particular heterogeneous samples were modeled by combining the contributions of several structures, each containing either one (periodic structure) or different layer types (mixed-layer structure -MLS) randomly interstratified (R=0). [18][19][20] These different contributions should be seen as a simplified way to describe the actual hydration heterogeneity of the sample under investigation, with the different layer types not being distributed at random in the different crystallites. However, the coexistence of these contributions does not imply the actual presence of populations of particles in the sample, as their relative proportions may vary as a function of RH for example. 18 To account for the heterogeneous distribution of the different layer types within smectite crystallites, layers exhibiting the same hydration state that are present in the different MLSs must have identical properties as they may be accounted for in one or the other structure depending on the RH. In particular for a given XRD pattern, each layer type must possess a constant crystal-chemistry in the different MLSs. It was possible to reproduce the profile of all experimental 00ℓ reflections of the experimental XRD pattern recorded on Ca-SWy-1 at 80% RH by considering two MLSs (Figure 1c) and the 2WS configuration for interlayer water. 18 Specifically, the position of the 005 reflection, the lowangle shoulder of the 002 reflection and the "tails" of the 001 reflection are satisfactorily reproduced by taking hydration heterogeneity into account. Accounting for hydration heterogeneity also helps reproducing the relative intensity of higher-angle reflections (002, 003, 004, and 005 for example) but significant discrepancies that could result from an incorrect structure model for interlayer water are still visible for high-angle reflections (Figure 1c). Specifically, the 006, 007 and 008 reflections are not satisfactorily reproduced, as for example the intensity ratio between the 007 and 008 reflections measured on experimental and calculated patterns are inconsistent. These discrepancies are reduced by increasing the Debye-Waller factor of H 2 O molecules (B wat ) from 2 to 11 Å 2 for this 2WS configuration of interlayer H 2 O molecules (Figure 1d). 18 However, such high values of the Debye-Waller factor are not sufficient to conceal the disagreement for the intensity ratio between 007 and 008 reflections, and thermal atomic fluctuations most likely do not adequately describe the positional distribution of H 2 O molecules in 2W smectite layers, and additional hypotheses have to be sought.

Interlayer configuration of 2W smectite layers as seen by Monte-Carlo

simulations. In the above calculations, H 2 O molecules are distributed in discrete planes, and the positional distribution of H 2 O molecules results only from their thermal motion. However, this simplified description of the smectite interlayer structure does not allow fitting the experimental XRD data (Figures 1c,1d) most likely because the description of H 2 O molecule positional disorder is incomplete. A more complete (realistic ?) description of the interlayer structure may be obtained from Monte-Carlo (MC) simulations which allow taking into account all interactions among interlayer species, as well as between these species and the 2:1 layer. 21 It is in particular possible to account for the hydration variability of interlayer cation which can form either inner-sphere or outer-sphere complexes with the 2:1 layer surface, leading to the existence or to the lack, respectively, of direct interactions with O atoms from the layer surface. In the latter case, these interactions are screened by H 2 O molecules from the cation hydration sphere. It has been shown that, as compared to other monovalent cations, K + cations tend to form inner-sphere complexes in montmorillonite interlayers and that these cations remain partially bound to the 2:1 clay surface even in the 2W state. 22,23 On the contrary, Li + and Na + cations in 2W smectites are located in the mid-plane of the interlayer. 23- The location of the layer charge deficit has also been shown to influence the hydration of interlayer Na + cations, the formation of inner-sphere complexes being favored by tetrahedral substitutions. 26 On the other hand, a majority of interlayer Na + cations is located in the midplane of the interlayer for octahedrally substituted 2W smectites. 23,24 A similar influence of the charge location was reported for K-and Li-saturated 2W smectites. 22,25 In contrast, whatever the charge location Mg 2+ cations are systematically octahedrally coordinated in 2W smectites and located in the mid-plane of the interlayer. 27,28 In any case, MC simulations most often indicate that H 2 O molecules do not form a discrete plane but rather show that they are distributed about a "most probable" position. In addition, the mixed charge location common in smectite layers, and more especially in those of natural samples, can lead to the coexistence in a single smectite interlayer of different complexes, thus broadening the water distribution profile by perturbing the hydrogen bond network and the orientation of the water dipole. 29 Even though MC simulations do not commonly account for smectite hydration heterogeneity, which is best revealed by XRD analysis, such a description of H 2 O molecules positional disorder could be the missing link toward a better structure determination of H 2 O configuration in 2W smectite layers.

Materials and Methods

Experimental. Samples investigated in the present work include two reference lowcharge montmorillonites (SWy-1 and SWy-2) available from the Source Clays Repository (http://www.agry.purdue.edu/cjohnston/sourceclays/index.html) and two synthetic saponite samples. The latter samples were selected because of their contrasting layer charges (0.8 and 1.4 per O 20 (OH) 4 ). 30,31 The size fractionation of all samples, and their homoionic saturation were performed as described by Ferrage et al. 18 For all samples, oriented slides were prepared by drying at room temperature a clay slurry pipetted onto a glass slide. XRD patterns were then recorded using a Bruker D5000 diffractometer equipped with a Kevex Si(Li) solid-state detector, an Ansyco rh-plus 2250 humidity control device coupled to an Anton Paar TTK450 chamber. Usual scanning parameters were 0.04°2θ as step size and 6s as counting time per step over the 2-50°2θ Cu Kα angular range. The divergence slit, the two Soller slits, the antiscatter and resolution slits were 0.5°, 2.3°, 2.3°, 0.5° and 0.06°, respectively. Data collection conditions (60 and 80% RH for Sr-saturated samples, 40 and 80% RH for Casaturated samples, and 80 or 90% RH for Na-saturated samples) were selected because of the high amount of 2W layers (>90%) present in these conditions.

18 Simulation of X-ray diffraction data. The algorithms developed initially by Drits and coworkers were used to fit experimental XRD profiles over the 2-50°2θ CuKα range using a trial-and-error approach. [START_REF] Drits | X-Ray structure analysis of mixed-layer minerals[END_REF][33][34] Instrumental and experimental factors such as horizontal and vertical beam divergences, goniometer radius, length and thickness of the oriented slides were measured and introduced without further adjustment. The mass absorption coefficient (µ * ) was set to 45 cm 2 g -1 , as recommended by Moore and Reynolds, 1 whereas the parameter characterizing the preferred orientation of the particles in the sample (σ*) was considered as a variable parameter. Additional variable parameters include the coherent scattering domain size (CSDS) along the c* axis which was characterized by a maximum CSDS value, set to 45 layers, and by a variable mean CSDS value (N). 35 In addition, because of the weak bonds between adjacent smectite layers, layer thickness was allowed to deviate from the average d 001 value. This cumulative deviation from periodicity, which is described as a "disorder of the second type", [START_REF] Guinier | Théorie et technique de la radiocristallographie[END_REF][START_REF] Drits | X-ray diffraction by disordered lamellar structures: Theory and applications to microdivided silicates and carbons[END_REF] is accounted for by introducing a variance parameter σ z . 18 z-coordinates of all atoms building up the 2:1 layer framework as well as those present in the interlayer of 0W and 1W layers were set as proposed by Moore and Reynolds. [START_REF] Moore | Jr X-ray Diffraction and the Identification and Analysis of Clay Minerals[END_REF] The interlayer structure of 2W layers has been refined to account for all features of experimental XRD patterns recorded on 2W-dominated samples. In particular, a double Gaussian distribution of H 2 O molecules along the c* axis (2WG) was assumed. This 2WG model accounts both for the presence of a unique plane of H 2 O molecules on either side of the mid-plane (Figures 1c,1d) and for the positional distribution of H 2 O molecules derived from MC simulations incomplete. The 2WG distributions considered in the present study are symmetrical relative to the interlayer midplane. They are characterized by the distance (∆d) between this mid-plane, where interlayer cations are supposed to be located, and the position of the maximum density of the Gaussian distribution. In addition, the total amount of interlayer H 2 O molecules was refined together with the full width at half maximum intensity (FWHM) parameter of the Gaussian distribution. In the resulting structure model, H 2 O molecules were introduced using a 0.05 Å step along the c* axis, with a B wat factor equal to zero, as thermal motion is taken into account in MC calculations.

Two parameters were used to assess the overall goodness of fit. The unweighted R p parameter was considered because this parameter is mainly influenced by the most intense diffraction maxima such as the 001 reflection which contains essential information on the proportions of the different layer types and on their respective layer thickness values. The R wp parameter was also used to better account for the overall fit quality, especially in the highangle regions. [START_REF] Howard | Modern Powder Diffraction[END_REF] Accessory quartz reflections were omitted for the calculation of these parameters. On their low-angle side, calculated XRD patterns are limited to ~5°2θ CuKα because significant discrepancies, possibly resulting from an incorrect description of crystalline defects not challenging the results described in the present study, 18 Additional details on the MC simulations can be found elsewhere. 25,[41][42][43] Density profiles determined from MC calculations for interlayer sodium and H 2 O molecules were introduced in the XRD profile calculation using a 0.075 Å step.

Results

Influence of the Gaussian distribution profile on the relative intensity of 00ℓ

reflections. Figure 2 illustrates the influence of the different parameters used to describe the assuming a periodic Ca-SWy-1 2W structure (layer thickness = 15.2 Å), and calculated intensities were systematically normalized to that of the 001 reflection. By increasing the total amount of H 2 O molecules the intensity of the 002, 003, 005 reflections greatly increases, that of the 008 reflection also increases but to a lower extent whereas 004, 006 and 007 reflections are essentially unaffected (Figure 2). As its influence on 007 and 008 reflection intensity is limited, the nH 2 O parameter will not affect significantly the intensity ratio between these two reflections which is a common and critical discrepancy between experimental and calculated profiles (Figures 1a-d). On the contrary, the 008:007 intensity ratio is strongly affected by the FWHM of the Gaussian distribution, this ratio being minimum for a Dirac distribution and increasing with the FWHM of the distribution. The 007 reflection is actually more intense than the 008 one for FWHM values larger than ~1.3 Å (Figure 2). In addition this parameter may be strongly constrained from its major influence on the intensity ratio between two intense reflections (003 and 005 reflections) which can be reversed by increasing the width of the Gaussian distribution of H 2 O molecules. However, the 003:005 ratio is also affected by the ∆d parameter which also affects the 008:007 intensity ratio, both ratios increasing with increasing ∆d values. By increasing either the ∆d parameter or the FWHM, the intensity of the 002 reflection is systematically decreased, whereas that of the 004 reflection is increased or decreased, respectively. The intensity calculated for the 006 reflection is low whatever the values used for these two parameters.

Modeling of XRD patterns.

For all XRD patterns recorded on smectite samples, calculations were performed using three different configurations of H 1d). The 2WS configuration provides a satisfactory fit to experimental patterns for 00ℓ reflections with ℓ<6. However, this model does not allow concealing the discrepancy observed over the high-angle range, and more especially for the 008:007 intensity ratio, even if the Debye-Waller factor of H 2 O molecules is maximized (B wat = 10-11 Å 2 ). [START_REF] Lipson | In International tables for X-ray crystallography[END_REF] In this case, the water content and the ∆d parameter are increased from 6.6 to 6.8 H 2 O per O 20 (OH) 4 and from 1.20 to 1.32 Å, respectively, as compared to the 2WS configuration with B wat = 2 Å 2 (Table 2). A 008:007 intensity ratio consistent with that observed experimentally can be obtained by considering the 2WG configuration for interlayer H 2 O molecules. In this case, broad Gaussian distributions were assumed (FWHM = 1.7 Å), and both the water content and the ∆d parameter were increased as compared to alternative interlayer configurations (Table 2). This 2WG configuration also allows fitting better the profile of the 005 reflection, but that of the 003 one is slightly altered as a result of a low-angle tail broadening (Figure 1e).

The combination of two structures, a main periodic one with only 2W layers and a second one containing the three layer types (Table 1), accounts for the hydration heterogeneity of sample Ca-SWy-2 at 40% RH, and leads to the coexistence of 2W, 1W, and 0W layers (95%, 4%, and 1%, respectively). 20 The 2WS configuration allows describing most features of the experimental XRD patterns (R p = 1.31% and R wp = 8.13% -Figure 3a).

However, the 008 reflection is significantly more intense than the 007 one. By increasing the Debye-Waller B wat factor from 2 to 11 Å 2 , the 008:007 intensity ratio appears closer to the experimental one although the two estimates of the fit quality are not affected (Figure 3b).

This ratio is best reproduced by assuming a 2WG distribution with a FWHM of 1.4 Å (Figure 3c) although R p and R wp parameters are almost unaffected. As compared to the 2WS mode, the total amount of H 2 O molecules in such 2WG configuration is considerably increased from 6.2 (assuming a B wat factor of 2 Å 2 ) to 7.8 per O 20 (OH) 4 (Table 2).

Sr-saturated montmorillonite. At both 60 and 80% RH, the hydration heterogeneity of sample Sr-SWy-1 is minimum as it contains an overwhelming proportion of 2W layers (95, and 96%, respectively -Table 1). 18 As for the Ca-saturated samples, the 2WS configuration for H 2 O molecules leads to a satisfactory fit to the experimental XRD patterns, especially for 00ℓ reflections with ℓ<6, and for the 008 reflection (Figures 4a,5a). However, significant discrepancies between experimental and calculated patterns are visible for the 002 reflection and the 008:007 intensity ratio. These discrepancies are significantly reduced by increasing the Debye-Waller factor of H 2 O molecules from 2 to 11 Å 2 , but they do not vanish completely (Figures 4b,5b). The optimum fit to the experimental XRD patterns was again obtained assuming a 2WG distribution of interlayer H 2 O molecules with a large FWHM value (1.2, and 1.5 Å for Sr-SWy-1 samples recorded at 60 and 80% RH, respectively -Table 2; Figures 4c,5c). For the two samples, both Rp and Rwp are lower for the 2WG configuration of interlayer H 2 O molecules than for the 2WS ones.

Na-saturated montmorillonite. At 80% RH, the Na-SWy-2 sample exhibits a high proportion (92%) of 2W layers whereas minor amounts of 1W and 0W layers (5%, and 3%, respectively) account for the hydration heterogeneity (Table 1). As for the previous sample, the 2WS configuration of H 2 O molecules leads to a satisfactory agreement between experimental and calculated data, especially for 00ℓ reflections with ℓ<6, and for the 008 reflection (Figure 6a). However, by using a Debye-Waller factor of 2 Å 2 the intensity calculated for the 006 and 007 reflections are too low as compared to the experimental one.

Increasing the B wat factor up to 11 Å 2 significantly reduces these discrepancies although the 008:007 intensity ratio remains imperfectly reproduced (Figure 6b). The optimum fit to the experimental data for this sample was again obtained assuming a 2WG distribution of H 2 O molecules (Figure 6c -R wp = 5.33%, R p = 2.59%). The ∆d and FWHM parameters of this 2WG distribution are 1.50 Å and 1.4 Å, respectively (Table 2). A similar fit to the experimental data (Figure 6d -R wp = 5.34%, R p = 2.57%) was obtained assuming the distribution of interlayer species shown on Figure 7, while all other parameters were kept constant (Table 1). This distribution was derived from the MC simulations performed using the NVT ensemble. MC calculated distributions exhibit a single peak for the oxygen atoms, and two for the hydrogen atoms, between the interlayer mid-plane and the surface of reported from IR spectroscopy results, 45 and from microscopic simulations. 24 As for all other samples, the amount of interlayer H 2 O molecules has to be increased, together with the ∆d parameter, as the positional distribution of these species increases (Table 2).

Na-saturated synthetic saponites. At 90% RH, the hydration heterogeneity of both synthetic saponites is minimum as they exhibit an overwhelming proportion of 2W layers (94, and 97% for Na-Sap 0.8 and Na-Sap 1.4 samples, respectively -Table 1). As compared to the natural ones, these two synthetic samples present larger CSDS along the c* axis, as evidenced by the sharpening of the 00ℓ reflections ( 8a,9a). Increasing the Debye-Waller factor up to 11 Å 2 leads to a perfect fit to the experimental data for the high-charge sample (Na-Sap 1.4 -Figure 9b),whereas significant discrepancies are still observed between experimental and calculated patterns for the low-charge sample (Na-Sap 0.8 -Figure 8b). For this latter sample, the optimum fit to the experimental data was again obtained assuming a 2WG distribution of H 2 O molecules in the smectite interlayer with ∆d and FWHM parameters (1.39 Å and 1.4 Å, respectively) similar to those obtained for natural samples (Figure 8c -Table 2). For the Na-Sap commonly used for XRD pattern simulations can be discarded as it systematically leads to major discrepancies between experimental and calculated profiles (Figure 10). [18][19][20] Specifically, the use of this usual configuration systematically leads to poor fits to the experimental XRD patterns for low-angle high-intensity reflections such as 003, 004, and 005 reflections (Figure 10). By contrast, the distribution of H 2 O molecules within a single plane on either side of the mid-plane interlayer (2WS configuration) allows both fitting the profiles and reproducing the relative intensities of the 00ℓ reflections with ℓ<6 (Figures 1c, 3a, 4a, 5a, 6a, 8a, 9a). 18 When assuming a Debye-Waller B wat factor of 2 Å 2 , this model leads to significant discrepancies for high-order 00ℓ reflections, which are partly resolved by 9c). In the 2WG model the interlayer cation is considered to lie in a fixed position located in the interlayer mid-plane and to have a Debye-Waller factor of 2 Å 2 . This hypothesis does not imply that the interlayer cations are not distributed as H 2 O molecules are, but it was assumed as a first approximation that thermal motion would be sufficient to account for their positional disorder In addition, the sensitivity to the positional disorder of these cations is much reduced as compared to H

2 O molecules as the former species accounts for a minor part of the overall electronic density in smectite interlayers. For example, at 80% RH Ca 2+ cations account for only 6% of the interlayer electrons (Table 3).

When comparing the electronic density due to interlayer H 2 O molecules deduced from MC calculations with that obtained from XRD profile fitting (Figure 11a), it is possible to note that the overall profiles are globally alike in spite of significant differences. In particular For the 2WG model, the contribution of interlayer H 2 O molecules to the structure factor of 00ℓ reflections for a periodic 2W smectite can be expressed as: 2) can be expressed as:

) 2 cos( )) 2 1 ( 2 cos( ) sin ( 4 ) 00 ( 00 2 ∑ ∆ ∆ - = m h z m m h d O H n f F l l l l π π λ θ (2)
) 2 cos( )) 2 1 ( 2 cos( ) sin ( 2 ) 00 ( 00 2 2 ∑ ∆ ∆ - = m h Z m m h d eff O H O H n f n F l l l l π π λ θ (3) where ) 2 cos( 2 2 ∑ ∆ = m h Z m m eff O H n n l π is the effective amount of interlayer H 2 O
molecules contributing to the structure factor. Equations ( 1) and ( 3 Thus both 2WS and 2WG models are essentially different although in both cases the contribution of interlayer H 2 O molecules to the structure factor is strongly decreasing with increasing ℓ indices. In the first case, the thermal motion of these interlayer species is entirely responsible for the decrease whereas in the latter model the decrease is related to the decreasing effective number of H 2 O molecules contributing to coherent diffraction effects.

Note that both models may produce similar diffraction effects if appropriate values are used for the parameters describing the positional disorder of interlayer molecules. However, unrealistically large values were obtained for the B wat parameter when fitting Na-Sap 0.8 (B wat = 30 Å2 ) as compared to Na-Sap 1.4 (B wat = 11 Å 2 ) although both XRD patterns were recorded under similar experimental conditions, and the 2WG configuration of H 2 O molecules appears as more realistic than the 2WS one. Additional support for the 2WG model arises from the close match between the number of interlayer H 2 O molecules determined using the 2WG model and that measured independently from water vapor isotherms. layer type by the relative abundance of this layer type and compared to that obtained from water vapor adsorption-desorption isotherm experiments (Table 3). 18 The water content determined by Ferrage et al. from XRD profile modeling assuming a 2WS model for the distribution of interlayer H 2 O molecules was reasonably consistent with that obtained from water vapor adsorption-desorption isotherm experiments. 14,16,18 However, the 2WG configuration provides the best agreement with the water contents determined experimentally from water vapor adsorption-desorption isotherm experiments, the XRD values lying most often between the values obtained on either branches of the isotherm (Tables 2,3).

FWHM of H 2 O molecule Gaussian distribution.

When using the 2WG model to describe the distribution of H 2 O molecules in 2W layers, the FWHM parameter represents the positional disorder of the species, which is characterized by the B wat factor in usual models.

One may note that the diffraction effects resulting from the two configurations are similar and lead to a significant decrease of the coherent scattering of H 2 O molecules with increasing diffraction angle (see above). However, the B wat factor should be about constant for a given species whereas the FWHM parameter can be structurally interpreted. For example, when increasing the RH, the FWHM of the Gaussian distribution systematically increases for Caand Sr-saturated montmorillonites (Table 2 -Figures 11b,11c) most likely to accommodate the steady addition of H 2 O molecules weakly bound to the interlayer cation. On the contrary, with increasing layer charge, Na-saturated saponite samples hold more H 2 O molecules for a given RH value in a narrower distribution (Figure 11d -Table 2). A possible origin for such narrowing of H 2 O molecule distributions is the increased polarization of these interlayer species resulting from a stronger undersaturation of surface oxygen atoms. 2). The ∆d values reported in the present study represent only indicative values that could be used for XRD profile modeling but a more complete study should be carried out to determine the key factors that influence this parameter.

Relative positions of interlayer cations and H

Consistency with reported interlayer structures of expandable 2:1 phyllosilicates.

Comparison with the present data. Among expandable 2:1 phyllosilicates, vermiculite and smectite are differentiated from their contrasting layer charge, vermiculite exhibiting a higher layer charge (1.2-1.8 per O 20 (OH) 4 ) than smectite (0.4-1.2 per O 20 (OH) 4 ). 46 This difference is usually revealed by the contrasting swelling behavior of the two minerals after magnesium saturation and glycerol solvation, vermiculite and smectite exhibiting basal spacings of ~14 Å and ~18 Å, respectively, after such treatment. [START_REF] Moore | Jr X-ray Diffraction and the Identification and Analysis of Clay Minerals[END_REF][START_REF] Calle | Hydrous Phyllosilicates (exclusive of micas)[END_REF]48 However, distinct hydration behavior has not been reported for these two mineral species, and the predominance of bi-hydrated layers has been documented for the two species as a function of relative humidity. As a consequence, these two expandable 2:1 phyllosilicates will be considered together in the following.

For modeling XRD results of clay minerals containing 2W layers, the interlayer water configuration usually assumed for bi-hydrated smectite is that used for the calculations showed on Figures 1a and 10 and already described (Type I -Figure 12). [START_REF] Moore | Jr X-ray Diffraction and the Identification and Analysis of Clay Minerals[END_REF] This model does not allow the description of experimental XRD patterns (Figure 10) and may be rejected.

Most of the three-dimensional structural determinations of 2W interlayer configuration were actually performed on vermiculite as this mineral frequently exhibits ordered stacking sequences and because its higher content of interlayer cations allows for a more accurate refinement of cation positions as compared to smectite. In addition vermiculite, as illite, presents an ordered distribution of interlayer cations which eases the structural characterization of the interlayer configuration as compared to smectite. 49 The structural studies devoted to the configuration of interlayer species have led to different structure models that will be described below.

In bi-hydrated Mg-vermiculite, Mg 2+ cations are located in the mid-plane of the interlayer with one sheet of H 2 O molecules on each side of this plane (Type II -Figure 12). 5,50,51,[START_REF] Bradley | Clays & Clay Minerals[END_REF][START_REF] Shirozu | [END_REF] According to this model, Mg is octahedrally coordinated by six H 2 O molecules whereas additional H 2 O molecules, which are weakly bound to the cation, are located on the same plane as the six cation-bound H 2 O molecules. [54][55][56][57] A Type II configuration of H 2 O molecules was also proposed for Na-saturated vermiculite, 58,59 and for Na-, Ca-and Li-rich altered phlogopites. 60 A second configuration of interlayer species has been proposed for Ca-saturated vermiculites (Type III -Figure 12). 58,61,62 In this model, two distinct coordinations are reported for Ca 2+ cations, two out of three Ca 2+ cations being octahedrally coordinated as in type II configuration, whereas remaining Ca 2+ cations exhibit a cubic coordination. This dual coordination induces the presence of two discrete planes of H 2 O molecules (planes 2 and 3 -Figure 12) in addition to that observed in the type II configuration, which holds most H 2 O molecules (plane 1 -Figure 12). The increased heterogeneity of H 2 O configuration in Ca-, Sr-, and Ba-saturated samples as compared to Mg-saturated ones was confirmed both from diffraction and IR results. 56,57 A Type III configuration of H 2 O molecules was also proposed for Na-saturated vermiculite. 58 Figure 13 compares the 2WG configuration of interlayer H 2 O molecules determined for Ca-SWy-2 (40% RH) in the present study with that reported in the literature for Ca-saturated vermiculite. 58,61 After normalization of the three distributions to the denser plane of H 2 O molecules, the three planes of H 2 O molecules appear closely related to the 2WG configuration proposed in the present study to describe the positional distribution of interlayer species.

To compare the ∆d values obtained in the present study with those reported in the literature (1.14-1.45 Å -Table 4), these values can be normalized to the thickness of the interlayer space to account better for the balance of the interactions with the interlayer cation on the one hand and the 2:1 layer on the other hand (Table 5). Following such a normalization procedure, the ∆d values determined for the 2WG configuration of H 2 O molecules are consistent with those reported in the literature whereas lower values are obtained when assuming a 2WS configuration.

In addition, z-coordinates along c* axis were recalculated together with typical distances between the 2:1 layer and the planes of H 2 O molecules, and between H 2 O molecules and interlayer cations (Table 4). For type II and III configurations the distance between the 2:1 layer and the densest plane of H 2 O molecules scatters between 2.36 Å and 2.82 Å and is consistent with the formation of H-bonds between interlayer H 2 O molecules and the clay framework. The distance between the densest plane of H 2 O molecules and the interlayer cation ranges from 1.14-1.45 Å.

Specific interlayer structure resulting from the presence of tetrahedral substitutions. A third configuration of water in 2W smectite has been envisaged for Na-beidellite samples, with Na + cations being partly engaged in the ditrigonal cavities of the 2:1 layers and the coordinated H 2 O molecules distributed on either side of the interlayer mid-plane which is devoid of atoms (Type IV -Figure 12). [START_REF] Ben Brahim | In 5th Meeting of the European Clay Groups[END_REF][START_REF] Ben Brahim | [END_REF]63 Such a migration of the interlayer cation from the interlayer mid-plane toward the 2:1 clay framework is consistent with MC simulations and IR spectroscopy results which both support the formation of inner-sphere complexes for monovalent cations in tetrahedrally substituted 2:1 phyllosilicates. 21,24,26,31 In the present study, similar distributions of interlayer species have been determined whatever the location of the layer charge deficit in agreement with previous reports of Type II and Type III configurations of interlayer species in tetrahedrally substituted 2W vermiculites. [58][59][60] The central location of Na + cations was found to be consistent with experimental XRD data even when Na + cations account for a significant part of the interlayer electronic density (13% of the interlayer electrons for sample Na-Sap 1.4 ). Furthermore, if a Type IV configuration is assumed for the distribution of interlayer species, significant discrepancies arise between experimental and calculated patterns, especially for the 002 and 003 reflections which are extremely sensitive to the presence of interlayer species at the interlayer mid-plane position (Figure 14).

In conclusion, the present data does not provide experimental evidence for the migration of monovalent cations toward the surface of tetrahedrally substituted 2:1 layers. 18 (c) Calculation performed accounting for hydration heterogeneities and assuming a 2WS configuration with B wat = 2 Å 2 for H 2 O molecules. 18 Hydration heterogeneity was described by assuming the coexistence of a major MLS containing 2W and 1W layers (95:5 ratio) and of a second structure containing the three layer types (2W:1W:0W = 85:13:2) in a 61:39 ratio (Table 2). ). e Mean thickness of the coherent scattering domain size along the c* axis (in layers). f Sigmastar parameter characterizing the sample orientation (in °). 1 g Standard deviation of the layer thickness parameter (in Å). 18 molecules distributed as one plane on either side of the interlayer mid-plane. The Debye-Waller temperature factor for water (B wat ) given in Å 2 . c 2WG corresponds to an interlayer configuration of H 2 O molecules distributed according to a Gaussian function on either side of the interlayer mid-plane. In this case, B wat = 0 Å 2 . d Water amounts determined experimentally from water vapor adsorption/desorption isotherms. Data are taken from Cases et al., 16 and from Bérend et al. 14 for divalent and monovalent cations, respectively. e The water contents are given in mmol of water per g of clay. f Personal communication from Laurent Michot (LEM, Nancy, France). a Configurations of interlayer water in bi-hydrated smectite layers are schematized on Figure 2. b Distances are measured in projection along the c* axis and given in Å. O layer , H 2 O, and cation stand for the outermost plane of oxygen from the 2:1 layer, the H 2 O molecules and the interlayer cations, respectively. c B wat is the Debye-Waller temperature factor reported for H 2 O molecules (in Å 2 ). d Debye-Waller factor is 11 Å 2 for plane (iii) and 2 Å 2 for planes (i) and (ii), respectively. e Basal distance d 001 along the c* axis is given in Å. f nH 2 0/nCat represents the ratio between the number of interlayer H 2 O molecules and that of interlayer cations. g Average value for the different planes of H 2 O molecules. h Distances are given for the denser plane of H 2 O molecules. a 2WS corresponds to an interlayer configuration of H 2 O molecules distributed as one plane on either side of the interlayer mid-plane. The Debye-Waller temperature factor for water (B wat ) given in Å 2 . b 2WG corresponds to an interlayer configuration of H 2 O molecules distributed according to a Gaussian function on either side of the interlayer mid-plane. c The distance along the c* axis between the interlayer mid-plane and the maximum density of the H 2 O molecule distribution (∆d) is normalized to the thickness of the interlayer (L. Tck. minus the thickness of the 2:1 layer -6.54 Å). d Average value calculated from the data reported for Type II and Type III configurations of interlayer species. 51,53,55,58-

Figure captions

  are observed over the low-angle region. 39 Monte-Carlo simulations. Monte-Carlo simulations in the NVT ensemble were used to obtain a detailed spatial distribution of the different species within smectite interlayers. The model montmorillonite-type smectite used in the simulations has a Na 0.75 (Si 8 )(Al 3.25 Mg 0.75 )O 20 (OH) 4 structural formula and exhibits substitutions only in the octahedral sheet. The simulation box includes two 2:1 layers, each consisting of 8 unit cells (total area: 20.72 Å × 17.94 Å, thickness of the 2:1 layer: 6.54 Å). The total negative charge of the 2:1 layers was thus compensated for by 6 Na + cations in the interlayer. The interlayer shift between adjacent 2:1 layers was set to different arbitrary values for the two interlayers considered and not allowed to vary during the calculation. For the typical layer thickness value (15.52 Å) determined for Na-montmorillonite by XRD profile modeling, the water content was estimated from the results of previous MC simulations performed with the NPT ensemble. Series of such simulations allows the determination of layer thickness as a function of water content, at constant pressure and temperature, 40 and the water content was found to be 9.5 H 2 O molecules per O 20 (OH) 4 . The resulting distributions of H 2 O molecules within 2W smectite interlayers were collected over 5 million MC steps, normalized and made symmetric with respect to the mid-plane of the interlayer. The 2:1 layers were considered as rigid, and modeled with the rigid SPC/E model (O-H bond 1.0 Å, angle H-O-H 109.47°, charges -0.848 e -and +0.424 e -for oxygen and hydrogen atoms, respectively). Applied interaction potentials were the Lennard-Jones 6-12 and Coulombic potentials. Each atom in the simulation cell was thus characterized by two Van der Waals parameters and by its charge.

  the 2:1 layer, and are characteristic of the presence of a single plane of H 2 O molecules on either side of the cation plane which is located in the center of the interlayer. The distance between the maximum of the oxygen distribution and the maximum of the hydrogen distribution closer to the 2:1 layer is ~1.0 Å which is the length of the O-H bond in the water molecule. This indicates a preferential orientation of the H 2 O molecules in the interlayer, with one of the O-H bonds almost perpendicular to the surface of the 2:1 layer. Similar configurations of H 2 O molecules in the interlayer of octahedrally-substituted smectites have been previously

  increasing the positional disorder of H 2 O molecules (B wat = 11 Å 2 -Figures 1d, 3b, 4b, 5b, 6b, 8b, 9b). However, except for sample Na-Sap 1.4 , such an increased B wat factor does not allow fitting satisfactorily the high-order 00ℓ reflections, which would require unrealistically high B wat factor values. In addition, the contrasting B wat factors adjusted for the two synthetic saponite samples recorded under similar RH conditions plead for a different origin to the actual positional disorder of H 2 O molecules in smectite interlayers. Distribution of H 2 O molecules according to a double Gaussian function. The 2WG model can be considered as an improved version of the 2WS model in which the actual positional disorder of H 2 O molecules is better accounted for (Figures 1e, 3c, 4c, 5c, 6c, 8c,

  the two planes of H 2 O molecules on either side of the interlayer mid-plane are much narrower in the MC calculations (FWHM ~ 0.7 Å as compared to ~1.4 Å for XRD profile fitting) which indicate also a significantly higher electron density in the interlayer mid-plane. The narrower distribution obtained from the MC simulation can be due in part to the fixed interlayer displacement between adjacent layers considered for the calculations although the influence of interlayer shift and/or layer rotation on the distribution of interlayer species derived from MC simulations is expected to be limited. The simple (simplistic ?) Gaussian functions used to model the distribution of H 2 O molecules are both shifted toward the interlayer mid-plane (by about 0.2 Å) and broadened as compared to MC calculations. Both the broadening and the shift of the Gaussian distributions are likely related to the specific profile of the MC distribution, and more especially to the high electron density in the interlayer mid-plane (Figure7). However, the XRD profiles calculated assuming the two models are almost identical (Figures6c, 6d), pleading for a limited sensitivity of calculated XRD patterns to these two parameters if the actual distribution profile is unknown.Validity of the 2WG configuration model. Similar XRD patterns may be calculated with 2WG and 2WS models by increasing the B wat factor in the latter model (Figures8c, 8d, and 9b, 9c). Because of the demonstrated sensitivity of calculated XRD patterns to the distribution of H 2 O molecules, this similarity can only result from similar contributions of H 2 O molecules to the structure factor in both models. Factors affecting the structure factor include the scattering power, the position and the amount of considered species. If the origin of the layer unit is set in the center of the layer octahedron, the contribution of H 2 O molecules to the structure factor of 00ℓ reflections for a periodic 2W smectite (2WS model) the layer thickness. ∆d is the distance between the interlayer mid-plane and the positions of the H 2 O molecules along the c* axis. With increasing values of ℓ, the contribution of H 2 of the thermal motion of H 2 O molecules.

  power of H 2 O molecules (B wat = 0), and ∆d is the distance along the c* axis between the interlayer mid-plane and the position of the maximum density of the Gaussian distribution. is the amount of H n m 2 O molecules at a given distance (m∆z, m being integer) from the maximum density of the Gaussian distribution. The sum equals the total number of interlayer H ∑ m n m positional distribution of H 2 O molecules disturbs their coherent scattering and thus decreases their absolute contribution to the structure factor. The decrease becomes more important as the ℓ indice increases. To quantify this decrease, Equation (

  2 O molecules. The distance (∆d) between the interlayer cations, which are located in the interlayer mid-plane, and the maximum density of the interlayer H 2 O molecule distribution function was also varied from one model to the other, the maximum ∆d values being obtained with the 2WG configuration of H 2 O molecules (Table

Figure 1 .

 1 Figure 1. Comparison between experimental and calculated XRD patterns for the Ca-

  (d) Calculation performed accounting for hydration heterogeneities and assuming a 2WS configuration with B wat = 11 Å 2 for H 2 O molecules. (e) Calculation performed accounting for hydration heterogeneities and assuming a 2WG configuration (see text for details).

Figure 2 .

 2 Figure 2. Relative intensities of 00ℓ reflections, after normalization to the 001 reflection, as a

Figure 3 .

 3 Figure 3. Comparison between experimental and calculated XRD patterns for the Ca-

Figure 4 .

 4 Figure 4. Comparison between experimental and calculated XRD patterns for the Sr-saturated

Figure 5 .

 5 Figure 5. Comparison between experimental and calculated XRD patterns for the Sr-saturated

Figure 6 .

 6 Figure 6. Comparison between experimental and calculated XRD patterns for the Na-

Figure 7 .

 7 Figure 7. Density profiles of interlayer species along the c* axis derived from MC

Figure 8 .

 8 Figure 8. Comparison between experimental and calculated XRD patterns for the Na-

Figure 9 .

 9 Figure 9. Comparison between experimental and calculated XRD patterns for the Na-

Figure 10 .

 10 Figure 10. Comparison between experimental XRD patterns and those calculated assuming the usual configuration of H 2 O molecules. 1 Hydration heterogeneity has been taken into

Figure 11 .

 11 Figure 11. Density profiles along the c* axis of the electron distribution in the interlayer of

Figure 12 .

 12 Figure 12. Schematic description of the different configurations proposed in the literature for

Figure 13 . 61 Figure 14 .

 136114 Figure 13. Comparison of the distributions of H 2 O molecules reported for bi-hydrated

  2 . b 2WG corresponds to an interlayer configuration of H 2 O molecules distributed according to a Gaussian function on either side of the interlayer mid-plane. c The number of H 2 O molecules is given per O 20 (OH) 4 . In this case, B wat = 0 Å 2 . d The distance, in projection along the c* axis, between the interlayer mid-plane and the maximum density of the distribution of H 2 O molecules (∆d) is given in Å. e The width of the Gaussian distribution of interlayer H 2 O molecules (FWHM) is given in Å.

Fig. 01 Ferrage

 01 et al. Fig. 02

Ferrage et al. Fig. 10 Ferrage

 10 et al. Fig. 11

  Ferrage et al. Fig.12

Fraction

  of interlayer thickness Normalized water content Ferrage et al. Fig. 13 Ferrage et al. Fig. 14

  O molecules distributed according to the 2WG configuration. Optimum parameters used to characterize smectite hydration heterogeneity, that is the relative proportions of the different MLSs coexisting in the sample and their compositions (relative proportions of 2W, 1W, and 0W layers) are reported in Table1together with the layer thickness values for the different layer types, N, σ*, σ z and the water content in 1W layers. For 2W layers, the water content, the ∆d parameter, and the FWHM of the Gaussian distribution are reported in Table2for the different configurations of interlayer H 2 O molecules.

2 O molecules in the interlayers of 2W layers: (i) a 2WS configuration with two planes of H 2 O molecules characterized by a B wat factor of 2 Å 2 and a ∆d parameter of 1.2 Å, 18 (ii) a similar 2WS configuration with a larger Debye-Waller factor (B wat = 11 Å 2 ), and (iii) a configuration with H 2 Ca-saturated montmorillonite. For sample Ca-SWy-1 at 80% RH, the calculations performed for 2WS configurations of H 2 O molecules (∆d = 1.2 Å) and B wat factors of 2 and 11 Å 2 have been described above (Figures 1c,

Table 1 -

 1 Figures 8, 9). Layer thickness of 2W layers decreases from 15.4 Å to 15.0 Å as the layer charge increases from 0.8 to 1.4 per O 20 (OH)4 (samples Na-Sap 0.8 and Na-Sap 1.4 , respectively -Table1). For both samples, the 2WS configuration of H 2 O molecules with B wat = 2 Å 2 allows fitting satisfactorily 00ℓ reflections with ℓ<6 (Figures

Discussion Shortcomings of the usual description of H 2 O molecule positional disorder in 2W smectite

  It should be noted first that the initial assumption of identical properties for all layers exhibiting the same hydration state and present in the different MLSs was verified for all samples, thus validating the proposed description of smectite hydration heterogeneity. In addition, the configuration of H

	1.4 sample, a fit
	similar to the one obtained with a 2WS distribution of H 2 O molecules and a high B wat factor
	was obtained assuming a 2WG distribution of H 2 O molecules (Figures 9b, 9c). However, the
	FWHM parameter of this distribution is significantly lower (0.8 Å) than those typically
	obtained for natural samples (1.2-1.7 Å -Table 2).

interlayers. By accounting for smectite hydration heterogeneity, it is possible to model experimental XRD patterns thus gaining additional insights into the structure of smectite interlayers. 2 O molecules within 2W smectite layers

Table 1 .

 1 Optimum structural parameters used for the simulation of experimental XRD profiles. Relative proportion of the different contributions to the diffracted intensity. b Relative proportion of the different layer types in the different contributions to the diffracted intensity. 2W, 1W, and 0W stand for bi-hydrated, mono-hydrated and de-hydrated smectite layers, respectively. c Layer thickness (L. Tck.) of the different layer types. d Number of H 2 O molecules in 1W layers (per O 20 (OH) 4

	Sample	Rel. ab. (%) a	2W b 1W b 0W b	L. Tck. 2W c	L. Tck. 1W c	L. Tck. 0W c	nH 2 0 1W d	N e	σ* f	σ z g
	Ca-SWy-2 (40%RH) h	87 100 0 13 60 30 10 0	15.18 12.60 10.00 3.2	8.7	6.5 0.35
	Ca-SWy-1 (80%RH) i	61 39	95 85 13 5	0 2	15.51 12.85 10.00 3.3	6.0	6.5 0.27
	Sr-SWy-1 (60%RH) i	82 100 0 18 75 15 10 0	15.53 12.58 10.00 3.5	7.5	5.5 0.35
	Sr-SWy-1 (80%RH) i	84 100 0 16 75 15 10 0	15.73 12.70 10.00 5.5	7.5	5.5 0.35
	Na-SWy-2 (80%RH)	90 10	96 60 30 10 2 2	15.52 12.55 9.60	3.2	8.2 11.0 0.25
	Na-Sap 0.8 (90%RH)	44 100 0 56 90 5	0 5	15.40 13.20 9.80	5.7 13.0 2.0 0.19
	Na-Sap 1.4 (90%RH)	91 100 0 9 70 20 10 0	15.00 12.90 9.80	5.0 12.0 11	0.12

a

  h Data from Ferrage et al.20 i Data from Ferrage et al.18

Table 2 .

 2 Structural parameters of the interlayer space determined from XRD profile modeling as a function of the assumed water configuration. 2WS corresponds to an interlayer configuration of H 2 O molecules distributed as one plane on either side of the interlayer mid-plane. The Debye-Waller temperature factor for water (B wat ) given in Å

	Sample	2WS, B wat = 2 a nH 2 0 c ∆d d	2WS, B wat = 11 a nH 2 0 ∆d	nH 2 0	2WG b ∆d FWHM e
	Ca-SWy-2 (40%RH)	6.2	1.20	6.6	1.30	7.8	1.34	1.4
	Ca-SWy-1 (80%RH)	6.6	1.20	6.8	1.32	10.0	1.37	1.7
	Sr-SWy-1 (60%RH)	6.0	1.20	6.8	1.32	8.5	1.40	1.2
	Sr-SWy-1 (80%RH)	6.0	1.20	7.0	1.41	9.5	1.52	1.5
	Na-SWy-2 (80%RH)	7.4	1.20	8.2	1.41	9.5	1.50	1.4
	Na-Sap 0.8 (90%RH)	8.5	1.20	9.3	1.33	10.5	1.39	1.4
	Na-Sap 1.4 (90%RH)	8.4	1.20	9.0	1.33	9.4	1.35	0.8

a

Table 3 .

 3 Optimum amounts of H 2 O molecules determined from XRD profile modeling for the different configurations of interlayer species, and from water vapor adsorption/desorption isotherms. Interlayer configuration of H 2 O molecules commonly used for the calculation XRD profiles including 2W layers. 1 b 2WS corresponds to an interlayer configuration of H 2 O

	Sample	Type I configuration a	2WS B wat = 2 b	2WS B wat = 11 b	2WG c	Ads./Des. d
	Ca-SWy-2 (40%RH)	7.25 e	8.32	8.58	10.11	8.62/10.17
	Ca-SWy-1 (80%RH)	7.06	8.30	8.54	12.36	12.85/13.70
	Sr-SWy-1 (60%RH)	7.06	7.60	8.85	10.72	7.87/9.29 (10.70/11.90)
	Sr-SWy-1 (80%RH)	7.14	7.69	8.95	12.09	9.83/10.45 (12.80/13.70)
	Na-SWy-2 (80%RH)	7.04	9.27	10.24	11.82	10.50/13.10
	Na-Sap 0.8 (90%RH)	7.00	10.62	11.60	13.07	13.39/14.49 f
	Na-Sap 1.4 (90%RH)	6.96	10.45	11.19	11.62	14.23/17.18 f
	a

Table 4 .

 4 Structural parameters of the different configurations reported in the literature for interlayer water in bi-hydrated smectite layers.

			Type I configuration a
	Reference	Sample	Cation-H 2 O (iii) b	Cation-H 2 O (ii) b	Cation-H 2 O (i) b	B wat d	c,
	Moore and Reynolds 1	2W-Smectite	1.20			1.06			0.35	11/2 §
			Type II configuration
	Reference	sample	O layer -H 2 O b	Cation-H 2 O b	d 001 e		nH 2 O/nCat f B wat	c
	Mathieson et al. 51 Mg-Vermiculite	2.76	1.14		14.34			-	5.4
	Shirozu et al. 53	Mg-Vermiculite 2.67 g	1.17 g		14.33			7.44	6.1
	Alcover et al. 55 Mg-Vermiculite	2.69	1.19		14.36			-	-
		Altered Ca-Phlogopite	2.77	1.41		14.96			8.60	-
	Le Renard et al. 60	Altered Na-Phlogopite	2.71	1.43		14.87		10.70	-
		Altered Li-Phlogopite	2.71	1.30		14.62			8.79	-
	Beyer et al. 59	Na-Vermiculite	2.70 g	1.44 g		14.85			4.00	3.9
			Type III configuration
	Reference	Sample	O layer -H 2 O b,h Cation-H 2 O b,h		d 001 e	nH 2 O/nCat f B wat	c
	De la Calle et al. 61	Ca-Vermiculite	2.78		1.45 f		14.92	7.34	5.5
	Slade et al. 58	Ca-Vermiculite Na-Vermiculite	2.82 2.66		1.41 f 1.42		14.89 14.85	8.02 5.58	2.5 3.9
			Type IV configuration
	Reference	Sample	O layer -H 2 O b	O layer -Cation b	Cation-H 2 O b		d 001 e	nH 2 O/nCat f B wat	c
	Ben Brahim et al. 13	Na-Beidellite	3.00	1.00		2.00	15.25	11.87	5

Table 5 .

 5 Distances along the c* axis between the interlayer mid-plane and the maximum density of the H 2 O molecule distribution normalized to the thickness of the interlayer space for the different configurations of interlayer species.

	Sample	2WS, B wat = 2 a	2WS, B wat = 11 a	2WG b
	Ca-SWy-2 (40%RH)	27.8% c	30.1%	31.0%
	Ca-SWy-1 (80%RH)	26.8%	29.4%	30.5%
	Sr-SWy-1 (60%RH)	26.7%	29.4%	31.1%
	Sr-SWy-1 (80%RH)	26.1%	30.7%	33.1%
	Na-SWy-2 (80%RH)	26.7%	31.4%	33.4%
	Na-Sap 0.8 (90%RH)	27.1%	30.0%	31.4%
	Na-Sap 1.4 (90%RH)	28.4%	31.4%	31.9%
	Mean value	27.1% ± 0.7%	30.3% ± 0.8%	31.8% ± 1.0%
	Literature mean value d		32.7% ± 2.1%

O molecules. For a given ℓ value, the

O molecules can be approximated by weighing the water content hypothesized for each

Acknowledgments

The results presented are a part of a Ph.D. thesis granted by Andra (French National Agency for Nuclear Waste Disposal). Andra is thanked for the permission to publish this manuscript and for financial support. BL acknowledges financial support from the CNRS/PICS709 program, and from the CNRS/SdU "postes rouges" fellowships granted to BAS. VAD and BAS are grateful to the Russian Science Foundation for partial financial support. Laurent Michot (LEM, Nancy -France) is thanked for the fruitful discussions about smectite hydration. Jean-Louis Robert (IST Orléans, France) kindly provided the synthetic saponite samples.

Literature cited


domain: info stat stat.ml

text:
Danya Bachir
email: danya.bachir@gmail.com

Ghazaleh Khodabandelou

Vincent Gauthier

Mounim El

Eric Vachon

Combining Bayesian Inference and Clustering for Transport Mode Detection from Sparse and Noisy Geolocation Data

Keywords: Mobile Phone Geolocation, Call Detail Records, Trajectory Mining, Transport Mode, Clustering, Bayesian Inference, Big Data

Large-scale and real-time transport mode detection is an open challenge for smart transport research. Although massive mobility data is collected from smartphones, mining mobile network geolocation is non-trivial as it is a sparse, coarse and noisy data for which real transport labels are unknown. In this study, we process billions of Call Detail Records from the Greater Paris and present the first method for transport mode detection of any traveling device. Cellphones trajectories, which are anonymized and aggregated, are constructed as sequences of visited locations, called sectors. Clustering and Bayesian inference are combined to estimate transport probabilities for each trajectory. First, we apply clustering on sectors. Features are constructed using spatial information from mobile networks and transport networks. Then, we extract a subset of 15% sectors, having road and rail labels (e.g., train stations), while remaining sectors are multi-modal. The proportion of labels per cluster is used to calculate transport probabilities given each visited sector. Thus, with Bayesian inference, each record updates the transport probability of the trajectory, without requiring the exact itinerary. For validation, we use the travel survey to compare daily average trips per user. With Pearson correlations reaching 0.96 for road and rail trips, the model appears performant and robust to noise and sparsity.

Introduction

The growing use of smartphones generates massive ubiquitous mobility data. With unprecedented penetration rates, mobile networks are supplying the largest geolocation databases. Mobile phone providers collect real-time Call Detail Records (CDR) from calls, text messages or data at no extra-cost for billing purposes. Still, traditional transport planning models have so far relied on expensive travel surveys, conducted once a decade. Consequently, surveys are rapidly outdated, while suffering from sampling bias and biased users' responses. Past research used CDR to estimate travel demand [START_REF] Wang | Estimating dynamic origin-destination data and travel demand using cell phone network data[END_REF], optimal locations for new transport infrastructures [START_REF] Berlingerio | Allaboard: a system for exploring urban mobility and optimizing public transport using cellphone data[END_REF], weekly travel patterns [START_REF] Calabrese | Estimating origin-destination flows using mobile phone location data[END_REF], activity-based patterns [START_REF] Jiang | Activity-based human mobility patterns inferred from mobile phone data: A case study of singapore[END_REF], urban land-use [START_REF] Jameson L Toole | Inferring land use from mobile phone activity[END_REF], impact of major events or incidents [START_REF] James P Bagrow | Collective response of human populations to large-scale emergencies[END_REF] and population dynamics [START_REF] Bachir | Using mobile phone data analysis for the estimation of daily urban dynamics[END_REF][START_REF] Khodabandelou | Population estimation from mobile network traffic metadata[END_REF]. A few studies used triangulation, based on signal strength e.g., in Boston U.S. [START_REF] Alexander | Origindestination trips by purpose and time of day inferred from mobile phone data[END_REF][START_REF] Wang | Transportation mode inference from anonymized and aggregated mobile phone call detail records[END_REF]. In Europe, privacy policies restrict triangulation usage to police demands. CDR and GPS data both respect privacy compliance for geolocation. Still GPS data collection requires users to install tracking applications and activate GPS, which has greedy battery consumption. Consequently, GPS samples represent subsets of users' trips while CDR generate locations from larger populations over longer time periods. However CDR geolocation is coarse, noisy and affected by the usage frequency of devices. Raw CDR provide approximate and partial knowledge of true users' paths, hence requiring careful pre-processing. Past methods on transport mode detection mainly involved GPS data and are hardly transposable to CDR. In addition, these studies applied supervised learning [START_REF] Gonzalez | Automating mode detection using neural networks and assisted gps data collected using gps-enabled mobile phones[END_REF][START_REF] Reddy | Using mobile phones to determine transportation modes[END_REF][START_REF] Zheng | Understanding transportation modes based on gps data for web applications[END_REF] requiring a training dataset of trajectories with transport mode labels. Transport modes were either collected via applications where users consent to enter their travel details, or manually identified using expert knowledge, which is a costly task. In real world scenarios, transport modes of traveling populations are unavailable. Therefore we need new unsupervised approaches to tackle this issue. This paper presents the first unsupervised learning method for transport mode detection from any CDR trajectory. As this is a first study, we focus on a bimodal separation between road and rail trips. In collaboration with a mobile phone provider, we process one month trajectories from the Greater Paris, which are anonymized and aggregated for privacy. Trajectories are represented as sequences of visited mobile network areas, called sectors. Our model combines clustering with Bayesian inference to determine the probability that cellphones traveled by road or rail knowing their trajectories on the mobile network. The transport probability of a trajectory is initialized with a prior obtained from the travel survey and updated with each new visited sector. Transport probabilities for sectors are derived after clustering sectors by transport type. Sectors features are constructed using both mobile networks and transport networks spatial properties. Then, for a subset of 15% sectors, we extract transport labels, being road or rail, (e.g., equipments inside train stations, on highways etc.) while the remaining sectors are multimodal. For each cluster, we use the binary labels to calculate continuous transport probabilities as the proportion of labeled sectors among total sectors. Trajectories are thus attributed the most probable mode among road, rail or mixed (i.e., when probabilities are close). For validation, we calculate daily average rail and road trip counts per user and obtain Pearson correlations with the travel survey above 0.96, for the 8 departments of the region. In the next sections, we review the literature in Sec. 2 and describe data engineering in Sec. 3. The methodology steps are presented in Sec. 4. Eventually, we discuss main results in Sec. 5 and provide conclusion.

Related Work

Common applications for geolocation data mining are the identification of travel patterns for personal travel recommendation [START_REF] Zheng | Learning transportation mode from raw gps data for geographic applications on the web[END_REF][START_REF] Zheng | Learning travel recommendations from user-generated gps traces[END_REF], anomalous behavior detection [START_REF] Xiaolin Pang | On detection of emerging anomalous traffic patterns using gps data[END_REF] and transport planning [START_REF] Jiang | Activity-based human mobility patterns inferred from mobile phone data: A case study of singapore[END_REF]. Several works used supervised transport mode learning from GPS trajectories. A multilayer perceptron was used to identify car, bus and walkers modes for 114 GPS trajectories in [START_REF] Gonzalez | Automating mode detection using neural networks and assisted gps data collected using gps-enabled mobile phones[END_REF]. Features were the average and maximum speed and acceleration, the total and average travel distance, the number of locations divided by travel distance and the number of locations divided by travel time. The best accuracy was 91% using a 10-folds cross validation. In [START_REF] Reddy | Using mobile phones to determine transportation modes[END_REF], speed and acceleration features were collected from 16 GPS trajectories. Several classification models (Decision Tree, Kmeans, Naïve Bayes, NNeighbor, SVM, Discrete and Continuous HMM) were compared. The Decision tree with Discrete Hidden Markov Model obtained the highest accuracy (74 %). Still, supervised approaches with GPS are constrained by the small size of the training data. Moreover, although transport labels can be collected for small GPS datasets, they are unavailable for CDR. Meanwhile, few studies tackled unsupervised transport mode detection. In [START_REF] Biljecki | Transportation mode-based segmentation and classification of movement trajectories[END_REF] fuzzy logic was used as a scoring function calculated between consecutive GPS traces. The transport score was calculated with boolean conditions on speed, distances to transport network and previous mode. Still, this work lacked a performance evaluation. In [START_REF] Nabavi Larijani | Investigating the mobile phone data to estimate the origin destination flow and analysis; case study: Paris region[END_REF], base stations located inside Paris underground were used to identify underground mode from CDR trips. A record detected by an underground antenna was labeled accordingly. This approach is limited as it relies exclusively on indoor equipment inside the underground. No additional modes were identified. To our knowledge, only one work addressed unsupervised transport mode learning for two modes, road and public transport, using triangulated CDR [START_REF] Wang | Transportation mode inference from anonymized and aggregated mobile phone call detail records[END_REF]. The approach applies travel times clustering followed by a comparison with Google travel times. Still, CDR low frequency induces important incertitude and delay on start and end travel times of CDR trips. Consequently a device may not be detected as traveling when the real trip begins and ends. Moreover the presented approach was demonstrated on one unique Origin and Destination (OD) pair which is not sufficient to validate the method. In dense urban areas, travel times can be affected by traffic states (e.g., rush hours) and can be identical for several modes, depending on the OD. Our work presents a novel method for transport mode detection by combining two unsupervised techniques, namely clustering and Bayesian inference. This model classifies millions of CDR trajectories into road and rail trips. Instead of clustering trajectories with features such as speed or travel time, highly impacted by the imprecision, sparsity and noise of CDR geolocation, we apply clustering on sectors and build spatial features using transport networks. A small subset of road and rail labels is collected for sectors in order to calculate sectors transport probabilities. After the Bayesian inference step, we conduct a large-scale validation for the complete region, using the travel survey. The high Pearson correlations, obtained on daily average trips per user, proves the method is generalizable, performant and robust to noise and sparsity.

Data Engineering

For this study, we collect anonymized CDR trajectories from the Greater Paris region, over one month. Sectors features are constructed using the base stations referential jointly with transport networks infrastructures. For data normalization, we introduce a specific procedure accounting for heterogeneous urban density. Label extraction is realized to gather transport labels for a small subset of sectors. For model validation we use the household travel survey from 2010 conducted by Île de France Mobilités-OMNIL-DRIEA [1]. Mobile providers do not have access to GPS coordinates of mobile phones. Although we know which base station is connected to a device, it is unlikely to encounter mobile users positioned exactly at the base station. Devices are located inside mobile network areas covered by base stations signal range. For this study, we use the mobile network referential of the Greater Paris region. This region has a 12000 km 2 area with more than 1200 cities and 12 millions inhabitants. It is covered by thousands of mobile network antennas. Each base station emits 2G, 3G or 4G radio signals. Cells are circular areas covered by signals (see Fig. 1). Each cell equipment is oriented toward one direction. The partitions of cells directions are called sectors. The average sector number per antenna is 3 where one sector covers 120 • around the base station. A cellular tessellation is composed of a multitude of overlapping areas. We use the sector tessellation to get rid of overlaps and create the voronoï partitions using sectors centroids (see Fig. 2). We associate each mobile phone record to a sector location.

Mobile Network

Transport Networks

Transport networks are used to construct sectors features. We retrieve rails infrastructures for underground, overground, tramway and train stations shared by Île-de-France Mobilité on their platform Open Data STIF [START_REF]Open Data STIF[END_REF]. In addition we use OpenStreetMap (OSM) [3] to retrieve highspeed rails and road infrastructures. Roads are categorized by traffic importance. We filter residential roads which have highest road count and lowest traffic.

Raw Features Construction

We construct our dataset D = {d rail , d road , d station , n road , n rail , w station } where features stand for:

d road : shortest distance between sector centroid and road (see Fig. 2).

d rail : shortest distance between sector centroid and rail network (see Fig. 2).

d station : shortest distance between sector centroid and train station centroid.

n road : number of roads intersecting the voronoi.

n rail : number of rail lines intersecting the voronoi.

w station : weight of train stations calculated as the sum of stations area intersecting the sector voronoi area.

Data Normalization

We aim to find transport mode usage in sectors. As our raw features are built with spatial information they are impacted by urban density. In the city center the density is higher than in the suburb. Consequently sector areas and distances to transport networks are smaller while there are more transport hubs. We normalize our features to reduce the bias induced by urban density over transport usage. We introduce a normalization specific to our problem:

d norm,m = d m ∑ i d i ∈ [0, 1] (1)
n norm,m = n m ∑ i n i ∈ [0, 1] (2)
w norm,station = w station A v ∈ [0, 1] (3)

Sector Label Extraction

A few base stations are located on transport hubs, such as rail lines, train stations, highways or tunnels. We process this information to construct labels for a small subset of antennas. We assume that each sector inherits from its base station label. We attribute rail labels to indoor equipments located inside the underground and train stations, which represent 4% sectors. We assign road mode to indoor antennas in tunnels, constituting less than 1% sectors. We add outdoor antennas on highways (11% sectors) to increase the size of the road subset. In total we obtain 15% transport labels. In what follows, we use our subset of sectors with categorical transport labels {road, rail}, as prior knowledge. Still, categorical transport labels are not appropriate for most sectors, including outdoor equipments. In urban areas, such as the Greater Paris, the classic scenario is to encounter several transport modes inside an outdoor sector because of mobile networks' coarse granularity. Thus, we aim to find continuous transport probabilities P ∈ [0, 1] for all sectors, where indoor labeled equipments have maximal probabilities P ∈ {0, 1}.

Trajectories Pre-processing

For this study, the mobile provider pre-processed raw anonymized users' positions using noise reduction and segmentation (see Fig. 3). For segmentation, users' locations were separated into stay points i.e., when users remain in the same area, and moving points i.e., when users are assumed traveling. We define a trajectory as a sequence of moving points T u j = {(X 0 , t 0 ), ..., (X l , t l )}, j being the j th trajectory of the user u. The i th position recorded at timestamp t i is X i = (x i , y i ), where (x i , y i ) are the centroid coordinates of the visited sector. One trajectory corresponds to one user trip. We construct 95 millions CDR trajectories from 2 millions anonymized users during one month. Similar trajectories are aggregated to respect privacy policies. In order to compare our results with household travel survey, which was conducted for residents of the Greater Paris region, the mobile provider filters users by home department (first two digits of billing address postcode) and exclude visitors.

Model

This section presents the unsupervised learning scheme combining clustering and Bayesian inference to estimate transport modes of CDR trajectories. First, the prior transport probability is obtained from the travel survey. Second, the transport likelihood is calculated from the observed records, such as each new visited sector updates the probability. In this perspective, we apply a clustering on sectors. Then, our subset of sectors labels is used to calculate transport probabilities within each cluster. Each sector is assigned a continuous score in [0, 1] reflecting the real transport usage inside i.e., the probability to detect more users on the roads or on the rails. For each trajectory, we assign the mode with highest probability. Eventually, results are validated against the survey.

Clustering

We aim to find transport clusters for mobile network sectors with an underlying hierarchical structure. Thus we use an agglomerative hierarchical clustering. The

Evaluation Metrics

We use internal evaluation metrics to assess the clustering performance and to identify the optimal cluster number. We used the Silhouette (S) to evaluate clusters separability [START_REF] Kaufman | Finding groups in data: an introduction to cluster analysis[END_REF] (see Eq. 4).

s ik = b(i) -a(i) max(b(i), a(i)) (4)
S k = 1 N k N k ∑ i=1 s ik (5)
S = 1 N ∑ k S k (6)
where a(i) is the average intra cluster distances for sector i and b(i) is the lowest value among average inter cluster distances. Here N k stands for the size of cluster k. The number of samples equals N. The optimal number of clusters K maximizes the silhouette [START_REF] Liu | Understanding of internal clustering validation measures[END_REF]. In addition we used the S dbw validity index.

S dbw (k) = Scat(k) + Dens db (k) (7)
where

Scat(k) = 1 k k ∑ i=1 σ(ν i ) σ(D) (8)
and

Dens db (k) = 1 k(k -1) k ∑ i,j=1 dens(u ij ) max(dens(v i ), dens(v j )) (9)
where ν i denotes centroid of cluster i and u ij is the middle point between clusters i and j i.e., at mid distance from the two centroids (ν i ,ν j ). The scattering index Scat is used to estimate the intra cluster compactness based on standard deviations σ of clusters over total dataset D. The term Dens db represents clusters densities. It calculates the average ratio of clusters middle point densities over clusters centers densities. The underlying assumption is that well defined clusters are denser around their centroids than at their mid distance. This index is a trade-off between clusters densities and variances. It has been depicted as the most performing among internal clustering evaluation metrics in [START_REF] Halkidi | Clustering validity assessment: Finding the optimal partitioning of a data set[END_REF][START_REF] Liu | Understanding of internal clustering validation measures[END_REF].

The optimal cluster number is found when the index reaches its minimum.

Probability Scores of Sectors Transport Mode

For each cluster k we calculate the score p k,m for transport mode m ∈ {rail, road}.

p k,m = N k,m N m ( 10
)
where N k,m is the number of labeled sectors belonging to class m in cluster k and N m is the total number of sectors from class m in the dataset. We normalize p k,m to obtain the probability P(m|S i ) ∈ [0, 1] of using mode m given a visited sector S i , belonging to a cluster k.

P(m|S i ) = p k,m ∑ j p k,j (11)
Unlabeled sectors obtain transport probabilities according to their cluster. In addition we update the probabilities of outdoor labeled sectors (i.e., highways) using Eq. 10 and 11. Indoor labeled sectors have binary probabilities in {0, 1}.

Bayesian Inference of Trajectories Transport Mode

Bayesian inference is used to determine the main transport mode associated to mobile phone trajectories. In this perspective, we calculate the probability

Clustering Evaluation

We first compare the three linkage types. Average and complete linkage fail to separate sectors in the city center, with any distance metric. One huge centered cluster is produced with tiny clusters located at the region borders. We retain ward linkage with euclidean distance which produce clusters of comparable size, evenly present across the region. In order to find the optimal number of cluster we draw the dendrogram of the ward agglomerative clustering (see Fig. 4). The latter shows k = 2 is a good cluster number as it corresponds to the highest distance gap between merges. A small k leads to a macroscopic partitioning. We look for a higher k to detect finer transport modes tendencies.

A clear cut was possible for k ∈ {3, 4, 5, 9}, which were therefore also good candidates. We decide to bound the cluster number between 2 and 10. We use additional intra-cluster metrics. We calculate S and S dbw with several k values (see Fig. 5). The silhouette reaches a maximum for k = 4, for which separability is the highest. According to the S dbw minimization criterion, the optimal number of clusters is k = 9, for which clusters are the most compact and dense. For k ∈ [5, 10] the silhouette a local maximum for k = 9. For our problem we favor the larger k hence we select k = 9. We visualize the 9 clusters with t-sne (see Fig. 6) and project them on the sectors map (see Fig. 7).

Sectors Probabilities and Visualization

We calculate the transport probabilities per cluster (see Tab. 1). We describe clusters regarding transport usage. Each cluster is displayed in Fig. 6 and Fig. 7.

Trajectories

We infer transport probabilities for one month trajectories, filtering bank holidays. We count the number of rail and road trips (see Fig. 8). Only 3% trips have probabilities in range [0.4, 0.6]. We consider such trips have mixed (or uncertain) mode. In Fig. 8 we observe hourly travel patterns for a week. For business days, peak hours occur in the morning and early evening, with a smaller midday peak at lunch time. Morning and evening peaks appear unbalanced. One reason is that mobile phone usage tends to be more important in the evening thus we detect more users and more trips. A second reason could be that users travel more at the end of the day. This phenomenon is more pronounced for road trips, the highest gap being on friday evening.

Comparison with Survey

We compare our results with the latest household travel survey, from 2010, for the Greater Paris. About 43000 residents were asked about their travels during their past day, outside holidays. We calculate mobility statistics from survey and MP results (see Tab. 2). We average survey trip counts per resident:

C S = ∑ k i=1 N i * w i ∑ k i=1 w i
where an individual i of weight w i reported N i trips for the day he was questioned. The weight w i was calculated during survey with socio-demographic information to rescale the individual to the entire population. Similarly we average CDR trip counts per day and per device:

C MP = ∑ U i=1 ∑ T t=1 1 U 1
T n u,i where U is the number of phones, T is the number of days and n i,t is the number of trips detected for phone i for day t. In the survey, transport modes are separated in two categories, motorized modes including public transport, cars and motorbikes, and unmotorized modes i.e., walk and bike. Our model outputs the majority mode of a given CDR trajectory, between rail and road. We first examine results for all residents (see Tab. 2). The survey indicates the average trip number per user during a business day is 4.16 for all modes and 2.45 for motorized trips. We found an average of 2.10 daily trips per person. It seems we were able to detect 86% motorized modes. Because of the coarseness of the mobile network, walkers might be considered as non moving as their movement occurs at a too microscopic scale. In addition, the detection of travels is affected by CDR frequency. When a device is turned-off or unused for a long period of time, users are undetected. Compared to the survey, 14% daily motorized trips are undetected in average. We further analyze results for residents aggregated by home given for the city center, first ring , second ring and department scale (first two digits of postcode). We calculate Pearson correlations between survey and CDR estimates for all trips, motorized, road and rail trips. In addition we calculate the ratio between road and rail trips:

C ratio = C road C rail
. There is a negative correlation between total survey trips and CDR trips, due to the possible undetection of unmotorized modes. Correlations for rail, road and ratio are all above 0.96 for the three rings scale and the department scale. Still we have smaller ratio than the survey. The department obtaining results most similar with the survey is the city center (Paris). For the latter we detect the same number of motorized trips. This means that all users' trips were detected, suggesting that mobile phone activity of travelers is more important in the city center. From these observations we emit several hypothesis to explain remaining differences. First, because of their cost, surveys are performed on small population samples. Despite the use of weights to scale the sample to the total population, results can still contain sampling bias in addition with users' responses bias. Second, travel surveys are performed every 10 years because of their high cost. The latest complete survey is anterior to our study (seven years difference) which can lead to differences in results. In particular, transport policies over were oriented to favor public transport in the Greater Paris (e.g., introduction of a unique price for transport pass that reduced the price for suburbs). This could have influenced users to take public transports, especially in the suburb. In our opinion trips segmentation might impact results. Indeed our trajectories are segmented based on stay times. Public transport users sometimes experiment waiting times in stations e.g., when users change lines, and signals loss when entering the underground. This could cause higher trip segmentation for CDR rail trips. At last we detect 100% trips in the city center versus 80% in the suburb. In parallel the city center has the highest rail transport usage. This could indicate a bias in mobile phone usage i.e., public transport users are more likely to call, text or navigate on the web than drivers. Therefore some road trips could possibly be undetected.

Conclusion

From mobile phone data mining we can capture travel behavior of urban populations on multimodal transport networks. Compared to traditional travel surveys, Call Detail Records are a low-cost and up-to-date knowledge base for smart transport research. In this paper, we have introduced a novel transport mode detection method using CDR trajectories from the Greater Paris. Our model uses three data sources: mobile network data, transport networks and household travel survey. After significant data pre-processing, we combine clustering on mobile network areas, called sectors, with Bayesian inference for trajectories. From the clustering we find 9 clusters best described transport usage in the region. Three clusters exhibit high road probabilities, two had high rail probabilities while four had mixed usage. We compare our final results on trajectories with the household travel survey. Trips are aggregated by users' home location, at the department scale. We calculate the average number of trips per day for each user, averaged over all users. We obtain Pearson correlations above 0.96 for motorized, rail and road modes. It seems we detect exclusively motorized trips, as walkers movements are too microscopic regarding the mobile network scale. To our knowledge this is the first method separating road from rail trips considering all CDR trajectories from all users, with substantial comparison with survey data. Still it is hard to obtain exact same results as the survey. First we might have a different trip segmentation.

When users travel, their path on the network are likely to be segmented into subtrips because CDR are affected by waiting times and signals loss. This phenomenon could be more pronounced for public transport travels, as users often change lines and wait in stations. In addition, the detection of travels is impacted by usage frequency of phones. We observe that trips are most likely to be undetected when road usage is predominant. At last, surveys might contain bias, be outdated and miss particular events. This makes validation a difficult task as no available data source is a perfect ground truth. Our work shows encouraging results yet we have several pending issues we want to address in future works. First, although our model proved to be robust to noisy locations, oscillations filtering could be enhanced during CDR pre-processing. Second, as our model outputs one dominant mode, we need to address multi-modal and uncertain behaviors. For future work, we will extend model evaluation with finer scale Origin-Destination trips. We look forward to adding a fourth data source (e.g., travel cards data) for validation. We aim to enrich our model with additional transport modes. Our final model will be implemented by the mobile phone provider for B-2-B with transport operators and urban planners.

Fig. 1 :

 1 Fig. 1: schema of a tri-sector antenna.The antenna is represented by the black dot. Circular areas are cells for 2G, 3G and 4G signals.

Fig. 2 :

 2 Fig. 2: Example of a voronoi sector and its associated shortest distance to transports axes. Five roads (colored lines) and one rail line (dashed line) intersect the sector.

Fig. 3 :

 3 Fig. 3: Transport mode detection workflow applied for this case study. Pre-processing steps annotated with (*) were performed by the mobile operator.

Fig. 4 :

 4 Fig. 4: Dendrogram for k ∈ [2, 10]. The xaxis is the height i.e., distances between clusters leaves and nodes. The yaxis shows the number of leaves per cluster.

Fig. 5 :Fig. 6 :

 56 Fig. 5: Silhouette (blue) and S dbw validity index (red) plotted in function of the number of cluster k

(a) C 1 9 Fig. 7 :

 197 Fig. 7: QGIS Clusters projection

  Fig. 7a. -C 3 , C 4 : rail dominated clusters with many underground sectors located in the city center. It corresponds to the red and yellow cluster on Fig. 7b. -C 5 , C 6 : mixed road clusters, shown in magenta and green on Fig. 7c. -C 7 , C 8 , C 9 : road clusters represented in black, orange and purple on Fig. 7d.

Fig. 8 :

 8 Fig. 8: Estimated trip counts are averaged per week day, per hour and per transport mode. Results are given for 1 month data from the Greater Paris.

  where d m ∈ {d road , d rail , d station }, n m ∈ {n road , n rail } and d norm,m , resp. n norm,m , is the normalized vector for feature d m , resp. n m . Feature w norm,station is the normalization of w station by voronoi area A v .

Table 1 :

 1 Transport Mode probabilities and cluster size for k = 9

	Cluster C 1	C 2	C 3	C 4	C 5	C 6	C 7	C 8	C 9
	Size (%) 14.7	8.50	12.4	4.67	2.20	10.5	24.4	5.60	17.1
	P RAIL 0.651	0.567	0.824	0.949	0.421	0.387	0.095	0.071	0.199
	P ROAD 0.348	0.432	0.176	0.051	0.579	0.613	0.905	0.929	0.801
	-C 1 , C 2 : mixed-rail clusters with a higher probability for rails, depicted in
	blue and cyan on

Acknowledgments

This research work has been carried out in the framework of IRT SystemX, Paris-Saclay, France, and therefore granted with public funds within the scope of the French Program "Investissements d'Avenir". This work has been conducted in collaboration with Bouygues Telecom Big Data Lab.

P(m|T u j ) to take a mode m ∈ {rail, road} knowing the trajectory T u j , using Bayes theorem: P(m|T u j ) = P(T u j |m) * P(m) P(T u j )

Trajectories are sequences of sectors {S 0 , ..., S l } visited by mobile phone holders. Thus we have P(T u j |m) = P(S 0 , ..., S l |m). We assume independence between sectors probabilities such as P(S i , S i+1 |m) = P(S i |m)P(S i+1 |m). This assumption is motivated by the need to reduce the computational cost of the calculation. Thus we can rewrite P(T u j |m) = ∏ l i=0 P(S i |m). Eq. 12 becomes:

The term P(m|S i ), previously calculated with Eq.11, is introduced by applying Bayes theorem a second time, to Eq. 12:

The term ∏ l i=0 P(S i ) P(T u j ) does not influence the mode choice. The prior transport probability P(m) can be seen as the initial guess, before observing records. The prior probability is obtained from the travel survey and is calculated as the average trip counts per user given the home location of cellphone holders, here at the department scale. For rail mode we have p rail,dep = AVG dep (c rail ) AVG dep (c rail )+AVG dep (c road ) ∈ [0, 1] and p rail,dep = 1p road,dep , where c rail and c road are the rail and road trip counts, for the day of survey, per user living in the department dep. At last we normalize the posterior transport probability to be in range [0, 1]. P(m|T u j ) ← P(m|T u j ) P(rail|T u j ) + P(road|T u j )

Finally we affect the mode obtaining the higher probability to each trajectory. When probabilities are in [0.4, 0.6] the mode is considered mixed.

Results

This section summarizes our main results. For the clustering we demonstrate how we determine the number of clusters. We describe clusters according to transport probabilities. From the Bayesian inference of trajectories' transport modes, we visualize transport flows per week day and observe the travel patterns. We provide detailed results comparison with survey, at department scale, using Pearson correlations as evaluation metric.