F Bouchelaghem
email: fatiha.bouchelaghem@upmc.fr

R Pusch

Fluid flow and effective conductivity calculations on numerical images of bentonite microstructure

Keywords: Hydraulic conductivity, Compacted water-saturated bentonite, Image analysis, Numerical Homogenization, Finite Element computations, Experimental comparison

Hydraulic conductivities of compacted water-saturated bentonite were computed based on the real microstructure. The Homogenization of Periodic Media approach employed fully acknowledges the heterogeneous and multiscale microstructure of clay, as well as locally varying physical flow properties. Consequently, three levels of description were considered : the microscopic level of clay particles, the mesoscopic level of clay aggregates, mineral grains and inter-aggregate porosity, and the macroscopic level of the sample subjected to fluid pressure gradients in the laboratory. Starting from the local description of fluid flow, the expression of the effective hydraulic conductivity tensor was derived. The soft and dense gels and the open voids may form a connected flow path or remain occluded. The local problems were solved on the microstructure obtained from a digitalized micrograph by image analysis. The contribution to macroscopic flow by the soft and dense gels was investigated in various configurations, and comparisons were made with hydraulic conductivity data for MX-80 bentonite.

Introduction

Clay microstructure

Smectites are characterized by a heterogeneous arrangement of nanometerscale lamellae stacked together to form clay particles, micrometer-scale grains, water and air which saturate the various pore spaces that appear on different scales [START_REF] Jozja | Etude de matériaux argileux albanais. Caracerisation multiéchelle d'une bentonite magnésienne[END_REF][START_REF] Holzer | 3Dmicrostructure analysis of hydrated bentonite with cryo-stabilized pore water[END_REF][START_REF] Keller | The pore structure of compacted and partly saturated MX-80 bentonite at different dry densities[END_REF][START_REF] Pusch | Microstructure of Smectite Clays and Engineering Performance[END_REF][START_REF] Tessier | Relation between the macroscopic behaviour of clays and their microstructural properties[END_REF]. As a result, clay microstructure controls most of their physical properties such as ion diffusion, swelling pressure and hydraulic conductivity. In water-saturated compacted bentonite, mesopores that contain free pore water may be distinguished from nanoscale pores that contain water molecules bound by surface forces. The multi-scale pore structure continuously changes during hydration, swelling and compaction processes, which complicates its characterization [START_REF] Holzer | 3Dmicrostructure analysis of hydrated bentonite with cryo-stabilized pore water[END_REF][START_REF] Tomioka | Analysis of microstructural images of dry and water-saturated compacted bentonite samples observed with X-ray micro CT[END_REF]. This explains the scarcity of microstructural data available to quantify the pore size distribution and pore connectivity in compacted bentonite. As pointed out in [START_REF] Holzer | 3Dmicrostructure analysis of hydrated bentonite with cryo-stabilized pore water[END_REF][START_REF] Keller | The pore structure of compacted and partly saturated MX-80 bentonite at different dry densities[END_REF], the proportions of free and bound porewater in compacted bentonite are the subject of ongoing discussions and investigations, and tomography techniques are still under development in order to quantify precisely the mesopore distribution in compacted swelling clays. In [START_REF] Pusch | Impact of microstructure on the hydraulic conductivity of undisturbed and artificially prepared smectitic clay[END_REF][START_REF] Pusch | The microstructure of MX-80 clay with respect to its bulk physical properties under different environmental conditions[END_REF] Transmission Electron Microscopy (TEM) micrographs of compacted and hydrated MX-80 bentonite displayed the formation of clay gels of variable density by linking of aggregates exfoliated from the expanding clay grains. Consequently, four main phases were distinguished: open voids which present no resistance to flow, clay gel-filled voids (designated as soft gels) which are characterized by a finite hydraulic conductivity, clay aggregates (designated as dense gels) and non-smectite grains which are assumed to be impermeable. Impervious zones and denser and softer gels were schematized by a grid of elementary cells with distinct hydraulic conductivities. The effective conductivity was then approximated by an analytical formula combining the individual cells in parallel or in series depending on their orientation with respect to the macroscopic flow direction. [START_REF] Tomioka | Analysis of microstructural images of dry and water-saturated compacted bentonite samples observed with X-ray micro CT[END_REF] used Xray microtomography for examining the morphological evolution of compacted Na-montmorillonite before and after saturation. They concluded that the outer montmorillonite sheets are likely to swell and form a gel that occupies the intergranular voids, whereas the inner sheets are not affected by the water saturation process. [START_REF] Holzer | 3Dmicrostructure analysis of hydrated bentonite with cryo-stabilized pore water[END_REF] and [START_REF] Keller | The pore structure of compacted and partly saturated MX-80 bentonite at different dry densities[END_REF] investigated the intergranular pore space (mesopores) of compacted and hydrated MX-80 bentonite using high resolution 3D imaging with Focused Ion Beam nanotomography. They observed that clay particles' hydration in non-compacted bentonite (dry density ρ dry = 0.39 g/cm 3 ) is associated with extensive exfoliation and dispersion of thin clay layers. At ρ dry = 1.24 -1.46 g/cm 3 , clay aggregates located in the intergranular pore space displayed a honeycomb structure of stacked clay platelets filled with a colloidal clay gel of lower density. At higher densities, the small remaining fraction of mesopores were elongated in shape and no longer inter-connected, which led the authors to conclude that fluid flow and transport cannot take place through the mesopores alone. In compacted hydrated bentonites such as those intended for radioactive waste disposal, the isolated mesopores are connected only through the nanopores present within the clay matrix. As we will see in Section 3.1 the macroscopic fluid flow in the case of unconnected mesopores is then controlled by the conductivity of the clay matrix.

Multiscale modelling of fluid flow in geomaterials

In recent studies [START_REF] Adler | Flow in simulated porous media[END_REF][START_REF] Al-Omari | Three-dimensional simulation of fluid flow in X-ray CT images of porous media[END_REF][START_REF] Blunt | Pore-scale imaging and modelling[END_REF][START_REF] Masad | Simulation of fluid flow and permeability in cohesionless soils[END_REF] the effective permeability of various porous materials (natural sandstones, carbonate reservoirs, asphalt pavements, sands, glass beads) was determined by solving the local equations of fluid flow. The microstructure was obtained directly from the sample by image processing, or idealized by utilizing statistical characteristics. The computed flow fields displayed preferential paths depending on the size and connectivity of the pores, with whole areas of the pores relatively isolated from the flow. In the study by [START_REF] Blair | Using two-point correlation functions to characterize microgeometry and estimate permeabilities of sandstone and porous glass[END_REF], a two-point spatial correlation function obtained from image processing was employed to estimate several microstructural parameters (porosity, specific surface area, mean grain size, effective pore size), and combined with a Kozeny-Carman relation to estimate the intrinsic permeability of porous glasses and natural sandstones. [START_REF] Blunt | Pore-scale imaging and modelling[END_REF] applied network modelling at larger scales to extract a topologically representative network consisting of pores and throats from the image, and the flow and transport properties were described semi-analytically. In [START_REF] Adler | Flow in simulated porous media[END_REF][START_REF] Al-Omari | Three-dimensional simulation of fluid flow in X-ray CT images of porous media[END_REF][START_REF] Masad | Simulation of fluid flow and permeability in cohesionless soils[END_REF], the computed velocity field was driven by the pressure difference prescribed between the inlet and the outlet of the microstructure. Periodic boundary conditions were imposed for the velocity components, Darcy's law was written at the macroscopic scale and the effective permeability tensor was simply defined as the ratio between the average fluid velocity vector and the known pressure gradient. The previous studies concerned materials of intrinsic permeability higher than 10 -12 m 2 , and did not investigate very low permeability materials such as compacted clays. Moreover, the pore space comprised only intergranular porosity with no nanoporosity such as encountered in clays. [START_REF] Ichikawa | Unified molecular dynamics and homogenization analysis for bentonite behavior : current results and future possibilities[END_REF][START_REF] Ichikawa | Seepage and consolidation of bentonite saturated with pure or saltwater by the method of unified molecular dynamics and homogenization analysis[END_REF]Ichikawa et al. ( , 2004a) ) and [START_REF] Wang | Numerical solutions for flow in porous media[END_REF] employed the Homogenization of Periodic Media to model fluid flow through clays. The microscopic and the macroscopic fluid flow equations were obtained from the Navier-Stokes equations using double scale asymptotic expansions, the two levels of description being related by the macroscopic permeability. Clay material was composed of lamellae of hydrated montmorillonite on the nanometric scale, while the mesoscopic domain consisted of clay particles and non-smectite grains. Water was distributed within the interlamellar and interparticle spaces, disregarding the mesopores. The computed effective permeability changed significantly with the void ratio, and comparison with experimental data [START_REF] Pusch | Waste Disposal in Rock[END_REF] confirmed the difficulties of assessing the permeability of very low permeability clays.

Modelling approach and outline of the paper

An experimental investigation on saturated bentonite [START_REF] Pusch | Impact of microstructure on the hydraulic conductivity of undisturbed and artificially prepared smectitic clay[END_REF] has been employed in a numerical homogenization model in order to relate the macroscopic hydraulic conductivity with structural and microtextural properties. With respect to existing works on compacted bentonites (Bouchelaghem and Jozja, 2009a,b;[START_REF] Ichikawa | Unified molecular dynamics and homogenization analysis for bentonite behavior : current results and future possibilities[END_REF][START_REF] Ichikawa | Seepage and consolidation of bentonite saturated with pure or saltwater by the method of unified molecular dynamics and homogenization analysis[END_REF]Ichikawa et al., , 2004a;;[START_REF] Wang | Numerical solutions for flow in porous media[END_REF] which are all based on simplified microstructures (Poiseuille-like flow for lamellae of hydrated montmorillonite, symmetrical arrangements of spherical or cylindrical grains on the mesoscopic scale), the macroscopic permeability tensor is computed on real microstructures, while both mesopores and non-smectite grains are taken into account, leading to new effective permeability formulations. According to [START_REF] Baltean | Etude de la dispersion d'un contaminant passif dans les milieux multiporeux[END_REF] and [START_REF] Wodie | Contribution a l'etude des milieux poreux par la methode de l'homogeneisation[END_REF], the effective permeability will be obtained on the macroscopic scale by upscaling from the mesoscopic scale alone, consisting here of the arrangement of clay aggregates, inter-aggregate voids and non-smectite grains. As mentioned in Section 1.1, microstructural data regarding the intergranular space [START_REF] Holzer | 3Dmicrostructure analysis of hydrated bentonite with cryo-stabilized pore water[END_REF][START_REF] Keller | The pore structure of compacted and partly saturated MX-80 bentonite at different dry densities[END_REF][START_REF] Pusch | The microstructure of MX-80 clay with respect to its bulk physical properties under different environmental conditions[END_REF][START_REF] Pusch | Microstructure of Smectite Clays and Engineering Performance[END_REF] allow us to distinguish between dense and soft parts of the heterogeneous clay matrix, implying regions of low permeability and highly permeable channels. Accordingly, we assume that the soft gels always contribute to the macroscopic flow, while the dense gels may be impervious to fluid flow, completely open to flow, or characterized by an effective conductivity that can be determined from a separate upscaling approach. The mesoscopic clay microstructure is described in Section 2.2, followed by the construction of a Finite Element mesh for numerical treatment in Section 2.3. Based on the description of effective fluid flow behaviour on the mesoscopic scale within clay aggregates (Section 2.4) and within inter-aggregate pores (Section 2.5), the method of asymptotic expansions is employed to derive the effective hydraulic conductivity tensor for unconnected mesopores (Section 3.1), for connected mesopores (Section 3.2), and when no mesopores are present (Section 3.3). Numerical simulations results and a comparison with experimental data on MX-80 bentonite are presented in Section 4. Finally, conclusions are drawn and future work is outlined in Section 5.

Formulation of the problem on the mesoscopic scale

Materials and sample preparation

MX-80 bentonite is a commercial smectite-rich clay, composed of numerous montomorillonite lamellae arranged in basic units or stacks [START_REF] Pusch | Microstructural evolution of buffers[END_REF][START_REF] Pusch | The microstructure of MX-80 clay with respect to its bulk physical properties under different environmental conditions[END_REF]. The chemical constitution is as follows: 65-75 % montmorillonite, 10-14 % quartz, 5-9 % feldspars, 2-4 % mica and chlorite, 3-5 % carbonates and chlorite, 1-3 % heavy minerals. The stacks are estimated to consist of 3-5 lamellae with equally oriented crystal axes in Na-saturated smectite, and about 10 lamellae if Ca is the dominant exchangeable cation [START_REF] Pusch | The microstructure of MX-80 clay with respect to its bulk physical properties under different environmental conditions[END_REF]. Buffer clay blocks are prepared by compacting air-dry MX-80 clay powder at a compaction pressure of 100 MPa. The dry density of the grains is ca 1980 kg m -3 at 10 % by weight and the dry density of the powder mass poured in the form and slightly compacted is ca 1200 kg m -3 . Transmission Electron Micrographs of very thin sections (about 3 10 -2 µm) are obtained by acrylate impregnation of the clay under confinement, the sample preparation is detailed in [START_REF] Pusch | The microstructure of MX-80 clay with respect to its bulk physical properties under different environmental conditions[END_REF].

Microstructure

The general description is based on observations of clay fabric that distinguish between two levels of arrangement: the particles and aggregation of particles [START_REF] Holzer | 3Dmicrostructure analysis of hydrated bentonite with cryo-stabilized pore water[END_REF][START_REF] Jozja | Etude de matériaux argileux albanais. Caracerisation multiéchelle d'une bentonite magnésienne[END_REF][START_REF] Pusch | Microstructure of Smectite Clays and Engineering Performance[END_REF][START_REF] Tessier | Relation between the macroscopic behaviour of clays and their microstructural properties[END_REF]. The microscopic domain of clay particles, which is characterized by the length l c , consists of stacked lamellae of montmorillonite, inter-layer water and lenticular pores (Bouchelaghem and Jozja, 2009a,b). Clay aggregates are represented by a periodical repetition of clay particles and inter-particle porosity. The mesoscopic domain is illustrated in Fig. 1a showing a TEM of MX-80 clay taken from [START_REF] Pusch | Microstructural evolution of buffers[END_REF]. The periodic microstructure Ω, of characteristic length l, is composed of the porous matrix Ω c of clay aggregates, water saturated inter-aggregate voids Ω f which may be connected or not, and non-smectite grains Ω s (in black in Fig. 1a). We will alternatively assume that the dense clay gels (in pink in Fig. 1a) are open to fluid flow and form a part of Ω f , or they will be considered as completely impervious, or they will be assimilated with a porous clay matrix Ω c of low hydraulic conductivity. In a similar way, the soft clay gels (in green in Fig. 1a) will either be considered as open to fluid flow (as part of Ω f ) or else be assimilated with a porous matrix of finite conductivity (contributing to Ω c ). The mesoporosity ϕ m is defined as the volume of inter-aggregate pores |Ω f | divided by the total volume |Ω|:

ϕ m = |Ω f | |Ω| = 1 - |Ω c + Ω s | |Ω| (1)
In the following, Γ cf and Γ cs respectively designate the interfaces between clay aggregates and inter-aggregate voids, and between aggregates and nonsmectite grains. Γ f s is the interface between mesopores and grains. Local periodicity is assumed for the microstructure and the physical properties, and Fig. 1b represents one quarter of the periodic cell Ω. The periodicity assumption is equivalent to the assumption of local statistical homogeneity and ergodicity for all the physical variables in random media. The following conditions of separation of scales are taken to be verified [START_REF] Bensoussan | Asymptotic Analysis for Periodic Structures[END_REF][START_REF] Sanchez-Palencia | Non Homogeneous Media and Vibration Theory[END_REF][START_REF] Rubinstein | Flow in random porous media: Mathematical formulation, variational principles, and rigorous bounds[END_REF][START_REF] Torquato | Random heterogeneous materials[END_REF]:

ε ′ = l c L ≪ 1, ε = l L ≪ 1, ε ′ ≪ ε (2)
where L designates the characteristic length of the macroscopic sample.

Image analysis and construction of the FE mesh

In order to create the geometry object to be meshed, the TEM micrograph in Fig. 1a is converted to a matrix format using the Image Processing Toolbox for Matlab (Matlab, 2010). In the resulting matrix of values ranging from 0 (black) to 255 (white), each level corresponds to a different density within the clay microstructure. A thresholding procedure based on a detailed exploration of the image is then performed, allowing to dissociate the three phases (soft gels and voids, dense gels, non-smectite grains) and resulting in a matrix with only three distinct values : 0 for the non-smectite grains, 125 for the dense gels and 255 for the soft gels and voids. The contours of the different phases are extracted and are imported into Comsol (Comsol, 2006) before being processed to obtain the surfaces occupied by the different phases. Fig. 1b represents the Finite Element unstructured mesh of the microstructure obtained by Delaunay triangulation. The relative fractions of the TEM picture are 41.6 % for soft gels and voids, 18.2 % for dense gels, and 40.2 % for non-smectite grains. The soft gels and voids occupy 44.6 % of the numerical microstructure, while the dense gels and grains represent respectively 20.2 and 35.2 % of the microstructure. The difference between the relative proportions of the phases in the micrograph and the FE mesh is mainly due to the noise which is present in the processed TEM picture and cannot be completely eliminated by the thresholding procedure. Another source of error is the presence of very small inclusions of one phase within another phase, that do not contribute to the overall fluid flow and are neglected during the mesh generation. By varying the digitalized microstructure resolution, Al-Omari and [START_REF] Al-Omari | Three-dimensional simulation of fluid flow in X-ray CT images of porous media[END_REF] and [START_REF] Blunt | Pore-scale imaging and modelling[END_REF] have observed that an approximately 1000-fold range from resolution to sample size may be necessary in order to obtain convergent permeability results. In our case, the high resolution employed, of the order of 10 nm/pixel in both the horizontal and the vertical directions, ensures a very good precision in capturing the actual microstructure and guarantees that the numerical permeabilities are convergent. Similar resolutions, comprised between 3.7 and 14.7 nm/pixel, have been employed to characterize the 2D and 3D mesopores distributions of MX-80 compacted hydrated bentonite [START_REF] Holzer | 3Dmicrostructure analysis of hydrated bentonite with cryo-stabilized pore water[END_REF]. Permeability computations similar to those depicted in Section 4 have also been performed on microstructures of increasing size, leading to convergent values for a microstructure size higher than two thirds of the microstructure illustrated in Fig. 1b.

Mesoscopic description of the domain occupied by clay aggregates

The derivation of the effective fluid flow behaviour through periodic media has been detailed in previous works for the case of constant fluid viscosity [START_REF] Baltean | Diffusion-convection in a porous medium with impervious inclusions at low flow rates[END_REF][START_REF] Baltean | Etude de la dispersion d'un contaminant passif dans les milieux multiporeux[END_REF][START_REF] Ichikawa | Unified molecular dynamics and homogenization analysis for bentonite behavior : current results and future possibilities[END_REF][START_REF] Ichikawa | Seepage and consolidation of bentonite saturated with pure or saltwater by the method of unified molecular dynamics and homogenization analysis[END_REF]Levy, 1987;[START_REF] Wodie | Contribution a l'etude des milieux poreux par la methode de l'homogeneisation[END_REF]. For random media, the homogenization procedure is similar to the homogenization of periodic media, provided that the assumption of periodicity is replaced by the assumption of local statistical homogeneity and ergodicity [START_REF] Rubinstein | Flow in random porous media: Mathematical formulation, variational principles, and rigorous bounds[END_REF][START_REF] Torquato | Random heterogeneous materials[END_REF]. The modelling approach proposed in Bouchelaghem and Jozja (2009a) in the case of variable viscosity can therefore be applied to random digitized clay microstructures such as the one illustrated in Fig. 1a.

In the inter-layer and inter-particle spaces Ω ′ f , fluid flow, which is assumed to be periodical, laminar and stationary, is governed by Navier-Stokes equations in which inertia effects are neglected. Fluid viscosity µ = µ(y ′ ) is taken to vary locally in the microstructure, in order to account for the "ice sheet" effect of viscosity increase near the particle surface [START_REF] Ichikawa | Unified molecular dynamics and homogenization analysis for bentonite behavior : current results and future possibilities[END_REF][START_REF] Ichikawa | Seepage and consolidation of bentonite saturated with pure or saltwater by the method of unified molecular dynamics and homogenization analysis[END_REF](Ichikawa et al., , 2004a;;[START_REF] Pusch | Microstructure of Smectite Clays and Engineering Performance[END_REF]. The mesoscopic description of flow in clay particles is obtained by performing an asymptotic analysis in terms of the parameter ε ′ [START_REF] Bensoussan | Asymptotic Analysis for Periodic Structures[END_REF]Levy, 1987;[START_REF] Sanchez-Palencia | Non Homogeneous Media and Vibration Theory[END_REF][START_REF] Rubinstein | Flow in random porous media: Mathematical formulation, variational principles, and rigorous bounds[END_REF], leading to (Bouchelaghem and Jozja, 2009a):

v 0 = -k j ∂p 0 ∂x j , p 1 = -p j ∂p 0 ∂x j (3)
where v 0 represents the first order approximation of fluid velocity, while p 0 = P (resp. p 1 ) is the first (resp. second) order approximation of fluid pressure, and k j = k j (x, y ′ ) verify the following system of equations:

µ 0 ∆ y ′ k j + ∇ y ′ µ 0 • (∇ y ′ k j + T ∇ y ′ k j ) + ∇ y ′ p j -e i = 0, ∇ y ′ • k j = 0 in Ω ′ f k j = 0 on Γ ′ sf (4)
where k j , p j are Ω ′ f -periodical. e i is the unit vector along axis y ′ i . The term

∇ y ′ µ 0 • (∇ y ′ k j + T ∇ y ′ k j ) expresses the contribution of variable viscosity.
By applying the averaging operator to Eq. ( 3)a, we obtain a generalization of Darcy's law at the first order of approximation:

v = 1 |Ω ′ | ∫ Ω ′ f vdV = (ε ′ ) 2 v 0 = -(ε ′ ) 2 k j ∂P ∂x j = -K∇ x P (5)
In Eq. ( 5), v = (ε ′ ) 2 v 0 represents the first order approximation of the fluid velocity within the porous space of clay particles, and K = (K ij ) is the hydraulic conductivity tensor of the porous clay matrix.

The problem description within the domain Ω c occupied by clay aggregates is completed by the continuity equation for the fluid phase, the non penetration condition on the interfaces Γ cs between clay aggregates and non-smectite grains (n is the outward unit normal on the interface), as well as by Ω-periodicity conditions:

∇ • v = 0 in Ω c (6) v • n = 0 on Γ cs (7) v, P Ω -periodical (8)

Mesoscopic description of the inter-aggregate voids

In inter-aggregate pores, under slow flow conditions of incompressible free water, transient and inertia terms are assumed negligible. According to nuclear magnetic resonance tests cited in [START_REF] Pusch | Microstructure of Smectite Clays and Engineering Performance[END_REF], the viscosity of water in pore spaces smaller than about 5 nm is higher than in free water. Some reduction in fluidity may also occur in voids as wide as 10 -2 µm, but for larger voids, the effect on porewater viscosity is estimated to be negligible. In the numerical simulations presented in Section 4, the smallest pore constrictions employed are around 0.05 µm. Consequently, fluid viscosity is assumed to be constant in mesopores and Navier-Stokes equations simplify to Stokes equations:

µ∆v f -∇p f = 0 in Ω f (9)
where µ, v f , and p f respectively designate fluid viscosity, velocity and pressure within the mesopores. Mass balance is written as follows:

∇ • v f = 0 in Ω f ( 10
)
Boundary conditions express fluid adherence on Γ f s and Ω-periodicity:

v f = 0 on Γ f s (11) v f , p f Ω -periodical (12)
We assume continuity of fluid pressures and normal velocities on both sides of Γ cf (n designates the unit outward normal on Γ cf ):

p f = P on Γ cf (13) v f • n = v • n on Γ cf (14)

Determination of the macroscopic hydraulic conductivity

Based on experimental evidence [START_REF] Holzer | 3Dmicrostructure analysis of hydrated bentonite with cryo-stabilized pore water[END_REF][START_REF] Keller | The pore structure of compacted and partly saturated MX-80 bentonite at different dry densities[END_REF][START_REF] Pusch | The microstructure of MX-80 clay with respect to its bulk physical properties under different environmental conditions[END_REF], the general description of a bentonite microstructure where clay gels of variable density partially fill the intergranular pore space has been retained. Three distinct formulations are obtained for the hydraulic conductivity at the macroscopic scale, depending on the distribution of the mesopores.

Effective hydraulic conductivity tensor with unconnected mesopores

Depending on the initial fabric of the bentonite, connected flow paths of inter-aggregate pores may not exist (Bouchelaghem and Jozja, 2009b;[START_REF] Souli | Physicochemical analysis of permeability changes in the presence of zinc[END_REF]. In that case, according to the local description given in Sections 2.4 and 2.5, fluid flow in the clay aggregates is described by Darcy's law Eq. ( 5) and the continuity equation Eq. ( 6). In water-saturated unconnected mesopores, fluid flow is governed by Stokes equations Eq. ( 9) and the incompressibility condition, Eq. ( 10). At the boundary between both domains, we assume continuity of fluid pressure Eq. ( 13) and normal fluid velocity Eq. ( 14) while, at the boundaries of solid grains, Eqs. ( 7) and ( 11) are valid. To complete the problem description, Eqs. ( 8) and ( 12) express Ω-periodicity. We must now consider that v and P depend on ε, in order to express the influence of mesopores and solid inclusions on fluid flow in the porous clay matrix.

As the unconnected voids Ω f are surrounded by the porous matrix Ω c , fluid flow in the mesopores is a consequence of fluid flow in the clay matrix, which implies that fluid velocities and pressures will be of the same order of magnitude in both domains. The asymptotic expansions are defined accordingly, each physical quantity being a function of both the macroscopic position vector x and the microscopic position vector y describing the local scale. The results of homogenization of the system of Eqs. ( 5) to ( 14), constituting a generalization of the results by [START_REF] Wodie | Contribution a l'etude des milieux poreux par la methode de l'homogeneisation[END_REF], are detailed by Bouchelaghem and Jozja (2009a). We find that the local velocity field in the porous clay matrix is given by the following expression at the first order of approximation:

v 0 i = -K 0 ij (δ jk + ∂Π k ∂y j ) ∂p 0 ∂x k (15)
Where the vector-valued function Π k does not depend on the macroscopic pressure gradient p 0 (defined everywhere in Ω), which enables us to formulate the equivalent behaviour under the form of a generalized Darcy's law. The macroscopic fluid velocity ⟨V 0 ⟩ is naturally defined as the volume average of the local fluid velocity fields in the clay aggregates Ω c and the mesopores Ω f :

⟨V 0 ⟩ = 1 |Ω| ( ∫ Ωc v 0 dV + ∫ Ω f v 0 f dV ) (16)
By use of the continuity equation, the adherence condition at the interface between the mesopores and the clay matrix, the periodicity condition and the divergence theorem, we obtain the final expression for the macroscopic velocity ⟨V⟩ at the first order of approximation. It defines the effective hydraulic conductivity tensor K ncv for the case of unconnected mesopores:

⟨V⟩ = (ε ′ ) 2 ⟨V 0 ⟩ = -K ncv ∇ x p 0 K ncv ij = 1 |Ω| ∫ Ωc K kl ∂ y k (y i + Π i ) ∂ y l (y j + Π j )dV (17)
In order to compute K ncv with Eq. ( 17)b from a volume averaging procedure, the following local problem for Π k (k = 1, 2) is solved by the Finite Element method (Comsol, 2006):

∂ ∂y i (K ij ∂Π k ∂y j ) = - ∂K ik ∂y i in Ω c Π k = -y k on Γ cf K ij ∂Π k ∂y j n i = -K ik n i on Γ cs (18) Π k = 0 on y k = 0, L e k ∂Π k ∂y j n j = 0 on y i = 0, L e k , i ̸ = k
where L e k (k = 1, 2) is the side length (along y k ) of the mesh in Fig. 1b.

Effective hydraulic conductivity tensor with connected mesopores

As in the case of unconnected mesopores, the mesoscopic description of fluid flow is given by the system of Eqs. ( 5) to ( 14). In connected mesopores, fluid velocity is now of order O(ε 2 ), while according to Eq. ( 5), fluid velocity in the domain of clay aggregates is of order O((ε ′ ) 2 ) and can be neglected. Consequently, the continuity of fluid velocities on Γ εcf will lead to a condition of null velocity for fluid flowing in the mesopores. The effective hydraulic conductivity tensor is obtained by upscaling of Stokes equations and neglecting fluid flow within clay aggregates, and the homogenization procedure leads to the local problem formulated in Section 2.4 in the particular case of constant viscosity:

v 0 f = -k f j ∂p 0 ∂x j ( 19
)
where p 0 is the macroscopic pressure, and the local problem is written as:

µ∆ y k f l + ∇ y p l -δ il e i = 0, ∇ y • k f l = 0 in Ω f k f l = 0 on Γ sf , k f l • n = 0 on Γ cf p l = 0, k f l • t = 0 on y l = 0, L e l ( 20
)
∇ y p l • n = 0, k f l • e j = 0 on y j = 0, L e j , j ̸ = l ∇ y k f l • n = 0 on y j = 0, L e j , j=1,2 k f l , p l being Ω ′ f -periodical.
e i is the unit vector along the local axis y i . Using Eq. ( 19) and defining the average fluid velocity ⟨V⟩ as the volume average of the local velocity field over Ω f , we have the following expression for ⟨V⟩ at the first order of approximation, defining the effective hydraulic conductivity K cv for connected mesopores as:

⟨V⟩ = ε 2 ⟨V 0 ⟩ = ε 2 |Ω| ∫ Ω f v 0 f dV = -K cv ∇ x p 0 , K cv ij = ε 2 |Ω| ∫ Ω f k f ij dV (21)
The local problem of Eqs. ( 20) to be solved for (k f j , p j ) is identical to the classical problem describing fluid flow of an incompressible Newtonian fluid under steady state conditions. The pressure gradient is taken equal to unity along the coordinate axis x l (of unit vector e l ), and zero in the other directions:

∂p 0 ∂x j = δ jl (22)
The effective hydraulic conductivity tensor K cv is then obtained through volumetric averaging on Ω f according to Eqs. ( 21) and ( 22).

Effective hydraulic conductivity with no mesopores

In highly compacted clays the mesoporosity is negligible (Bouchelaghem and Jozja, 2009b;[START_REF] Holzer | 3Dmicrostructure analysis of hydrated bentonite with cryo-stabilized pore water[END_REF][START_REF] Keller | The pore structure of compacted and partly saturated MX-80 bentonite at different dry densities[END_REF], and Eqs. ( 5) to ( 8) describe fluid flow in the clay aggregates. v and P depend on ε in order to account for the influence of solid inclusions at the mesoscopic scale on fluid flow in the porous clay matrix. This problem has been studied by [START_REF] Baltean | Diffusion-convection in a porous medium with impervious inclusions at low flow rates[END_REF]; [START_REF] Baltean | Etude de la dispersion d'un contaminant passif dans les milieux multiporeux[END_REF] and [START_REF] Wodie | Contribution a l'etude des milieux poreux par la methode de l'homogeneisation[END_REF]. By introducing asymptotic expansions, and using the same approach as described in Section 3.2, we obtain the local fluid velocity as:

v 0 i = -K ij (δ jk + ∂Π ′ k ∂y j ) ∂p 0 ∂x k (23)
The macroscopic fluid velocity is defined as the volume average of the local velocity field v 0 , and from Eq. ( 23) we obtain at the first order of approximation:

⟨V⟩ = (ε ′ ) 2 ⟨V 0 ⟩ = -K nv ∇p 0 , K nv ij = 1 |Ω| ∫ Ωc K ik ∂(y j + Π j ) ∂y k dV (24)
Eq. ( 24)b defines the effective hydraulic conductivity tensor K nv in the particular case when there are no mesopores present in the microstructure, and the system to be solved for Π ′ k takes the following form:

∂ ∂y i (K ij ∂Π k ∂y j ) = - ∂K ik ∂y i in Ω c K ij ∂Π k ∂y j n i = -K ik n i on Γ cs (25) Π k = 0 on y k = 0, L e k 2 ∂Π k ∂y j n j = 0 on y i = 0, L e i 2 , i ̸ = k
Knowing Π k , the effective hydraulic conductivity K nv is obtained by Eq. ( 24)b through volume averaging.

Results and discussion

Numerical computations have been performed in the case of isotropy of the clay matrix : K ij = K c δ ij , where K c is constant. The longitudinal permeability corresponds to the horizontal axis x while the transverse permeability is associated with the vertical axis y in the system of cartesian coordinates.

Connected mesopores

This situation is encountered for the determination of the transverse hydraulic conductivity, and either the connected mesopores correspond to the soft gels alone and the contribution of dense gels is neglected (Case a.), or the mesopores are associated with both the soft and dense gels (Case b.). Case a. Soft gels. The local problem of Eqs. ( 20) is solved in the domain occupied by the soft gels only, while the dense gels are assumed to be impermeable to flow. Fig. 2 illustrates the norm of the fluid velocity field computed within one quarter of the periodic cell. Fluid flow is highly localized and limited to the connected gel network at the left in the microstructure. The transverse intrinsic permeability value is equal to K cv yy = 1.041×10 -17 m 2 . The range of corresponding hydraulic conductivity values, obtained from the classical relationship K hydr yy = ρ bulk g µ K cv yy , with ρ bulk = 1300 kg/m 3 [START_REF] Pusch | Impact of microstructure on the hydraulic conductivity of undisturbed and artificially prepared smectitic clay[END_REF] and µ confined between 1.8 and 4 mPa s [START_REF] Paumier | Facteurs determinant la rheologie du systeme argile-eau pour des suspensions de smectites[END_REF], are listed in Table 1. In the horizontal direction there are no connected flow paths for the soft gels alone, therefore the effective longitudinal permeability cannot be computed using this approach, and will be obtained instead under the assumption of unconnected mesopores. The off-diagonal term K cv xy = K cv yx = 5.387 × 10 -18 m 2 , obtained from Eq. ( 21), is not negligible compared to the transverse permeability.

Case b. Soft and dense gels. When all the gels are assumed to be open to fluid flow, the computed transverse effective permeability is equal to K cv yy = 3.787 × 10 -16 m 2 . Using the same data as in Case a., the corresponding range of values for K hydr yy are listed in Table 1. The fluid velocity field is represented in Fig. 3, which shows connected flow paths as well. As in the simulations of [START_REF] Masad | Simulation of fluid flow and permeability in cohesionless soils[END_REF], the fluid flow is locally laminar but irregularly distributed within the void space. Narrow necks control the flow, while large portions of soft gels do not contribute to the overall fluid flow through the microstructure: fluid flows along preferential paths.

In the longitudinal direction there are no connected flow paths, the corresponding effective conductivity will be computed as below in Section 4.2. The off-diagonal component of the permeability tensor is obtained by integrating the horizontal velocity over the whole domain of soft and dense gels. From Eq. ( 21), we deduce K cv yx = 1.0923 × 10 -16 m 2 .

Unconnected mesopores

On the micrograph displayed in Fig. 1a, the soft gels alone do not form a continuous network throughout the whole microstructure. In that case, the hydraulic conductivity is obtained by solving the local system of Eqs. ( 19) in the domain Ω c occupied by the dense gels of hydraulic conductivity K c . Fig. 4 represents the distribution of fields Π 1 and Π 2 . Dirichlet boundary condition on the interface between the soft and dense gels (Π k = -y k on Γ cf , k = 1, 2) plays the role of a source term, (Bouchelaghem and Jozja, 2009a). The numerical intrinsic permeabilities are given in Table 1. We can also compute the off-diagonal components of the permeability tensor from Eq. ( 17), K ncv xy = 0.00197K c , and K ncv yx = 0.0020K c , showing, as expected, that the effective permeability tensor is symmetric. After diagonalization, the principal values of the permeability tensor are given by: K 1 = 0.1098Kc, K 2 = 0.0947Kc, these results are very close to the permeabilities computed in the cartesian coordinate system, the principal directions being only slightly inclined (about 0.131 rad) with respect to the cartesian axes.

No mesopores

The permeabilities are obtained by solving the system of Eqs. ( 26) and the averaging defined in Eq. ( 24). We distinguish two cases, depending on whether the dense gels contribute to the overall fluid flow or not. Case a. Soft gels. The transverse and longitudinal permeability values are summarized in Table 1. Fig. 5 illustrates the local distribution of fields Π 1 and Π 2 . As in the unconnected macrovoids case, Neumann condition (26)b on the interface between the clay gels and non-smectite grains ( ∂Π k ∂yi n i = -n k on Γ cs , k = 1, 2) plays the role of a linear source of force. Case b. Soft and dense gels. The computed transverse and longitudinal permeabilities are listed in Table 1. Fig. 6 displays the distributions of Π 1 and Π 2 within the microstructure. As in Case a., the longitudinal permeability is negligible with respect to the transverse permeability. The transverse permeability is of the same order of magnitude as the value obtained when soft gels are treated as unconnected mesopores.

Discussion

The permeability values computed for 3D microstructures are always higher than those obtained by a 2D analysis. According to [START_REF] Adler | Flow in simulated porous media[END_REF]; Al-Omari and [START_REF] Al-Omari | Three-dimensional simulation of fluid flow in X-ray CT images of porous media[END_REF] and [START_REF] Masad | Simulation of fluid flow and permeability in cohesionless soils[END_REF], the discrepancy between 2D and 3D computations can be attributed to the inherent limitation of a 2D computation in capturing the real 3D flow path in a highly heterogeneous material such as sand or clay, as well as to the difference in specific surface area between the discretized contours of the pore space and the real microstructure. Therefore, 2D computations can be improved by increasing the resolution of the images in order to capture the real specific surface. As the surface area of the sand samples investigated by [START_REF] Masad | Simulation of fluid flow and permeability in cohesionless soils[END_REF] increases, the number of flow paths decreases and they are more tortuous as well, while narrower necks appear within the porous space leading to a ratio between the average 3D permeability and the average 2D permeability of 1.28. Since in the present study there is no significant difference between the experimental and numerical surface areas, and as the surface area is notably higher for clays than for sands, we may reasonably assume that the 2D permeability is of the same order of magnitude as the 3D permeability.

Another aspect of interest is the anisotropy of the computed permeability values, which is measured by [START_REF] Masad | Simulation of fluid flow and permeability in cohesionless soils[END_REF] as the ratio r k of the horizontal permeability k h to the vertical permeability k v , r k = k h kv . In their study, r k ranges from 1.1 for spherical glass beads to 1.67 for silica sand, and the overall agreement between the 2D and 3D computations, as well as between the computations and the experimental values, increases with r k . [START_REF] Masad | Simulation of fluid flow and permeability in cohesionless soils[END_REF] have shown that the permeability anisotropy is related to an increase in surface area and to the anisotropy of the porosity distribution in the samples. Our results are consistent with their findings. The anisotropy ratio for the case of unconnected intergranular voids is r k = 0.8624. In the connected mesopores as well as in the no-mesopores simulations, the permeability is highly anisotropic since the longitudinal permeability is negligible, implying that r k ≃ 0.

Comparison with experimental data

MX-80 in Ca-form. For Wyoming clay in Ca-form, the gel conductivity value K gel = 7 ×10 -5 m/s has been obtained by [START_REF] Pusch | Microstructure of Smectite Clays and Engineering Performance[END_REF] for ρ bulk = 1570 kg/m 3 . Using this value for K c , we obtain in the unconnected mesopores simulations the results synthesized in Table 2. Results are also given for the simulations with no mesopores, with a permeable fraction limited either to soft gels alone, or attributed to both the soft and dense gels. No mesopores Soft gels 1.505 ×10 -6 ∼ 0

No mesopores Soft + dense gels 9.058 ×10 -6 3.4767 ×10 -9

Compared to the experimental bulk conductivity of K Ca exp = 2 ×10 -9 m/s, the model with connected mesopores and permeable soft and dense gels behaves well, and the no-mesopores model in the longitudinal direction gives interesting results as well. For the other models, the computed permeability is either several orders of magnitude higher than the experimental value or negligible. More values are required to draw a definite conclusion. Still, the calculations indicate that both the soft and dense clay gels contribute to flow for MX-80 in Ca-form. MX-80 in N a-form. For Wyoming clay in N a-form, the gel conductivity value K gel = 2 ×10 -10 m/s given in [START_REF] Pusch | Microstructure of Smectite Clays and Engineering Performance[END_REF] is four orders of magnitude lower than the Wyoming clay conductivity in N a-form measured under similar compaction (ρ bulk = 1570 kg/ m 3 ). Using this value for K c , the results of the simulations are gathered in Table 3. The experimental conductivity measured in [START_REF] Pusch | Microstructure of Smectite Clays and Engineering Performance[END_REF] is equal to K N a exp = 8 ×10 -11 m/s. The best agreement between the model and the experimental value is when the soft gels are assumed to make up the connected mesopores network. When both the soft and dense gels are included in the connected mesopores model, the computed permeability values are too high. The permeabilities computed using unconnected mesopores and no mesopores are slightly lower but of the same order of magnitude as the experimental data.

Using idealized microstructures, Bouchelaghem and Jozja (2009a) obtained an overall agreement between the experimental results of [START_REF] Pusch | Microstructure of Smectite Clays and Engineering Performance[END_REF] and the model with connected mesopores. However, the modelling with the mesopores filled with a permeable fraction of gels of finite conductivity must not be dismissed. Bouchelaghem and Jozja (2009b) have shown that for watersaturated Wyoming bentonite in Mg form, numerical conductivities using unconnected spherical mesopores in 3D microstructures compare very well with the conductivities measured before and after pollutant permeation.

Conclusion

In the present study, the effective hydraulic conductivity tensor was computed by solving the local Stokes equations of fluid flow through a 2D microstructure of a hydrated MX-80 bentonite. Starting from a digitized TEM micrograph that distinguishes between open voids, soft and dense gels and nonsmectite grains, different models were proposed which differ by the way the clay gels are treated. Those models have been developed in order to evaluate the contribution to the overall hydraulic conductivity of clay gels of varying density, and to assess the effect of mesopore connectivity. The different models were compared with existing experimental data available for MX-80 clay under Ca and N a form, showing that the description with connected mesopores does not necessarily lead to the highest values of hydraulic conductivity. More specifically, for clay under Ca-form, the clay gel hydraulic conductivity is much higher than the bulk hydraulic conductivity measured, implying that the computed values of hydraulic conductivity with the unconnected mesopores model are too high. In that case, the best agreement with the experimental conductivity was obtained in the connected mesopores configuration, when both the soft and dense gels are open to fluid flow. For clay under N a-form, the hydraulic conductivity computed with connected mesopores associated with soft gels alone was close to the experimental value. The latter results are consistent with the experimental observations by [START_REF] Pusch | The microstructure of MX-80 clay with respect to its bulk physical properties under different environmental conditions[END_REF] showing that in the Ca-state the smectite stacks are thicker, the soft clay gels consist of bigger aggregates yielding larger open voids than in the N a-state. However, for clay under N a-form, the models with unconnected mesopores or with no mesopores resulted also in numerical hydraulic conductivities of the same order of magnitude as the experimental hydraulic conductivity.

A constant value was employed for the hydraulic conductivity of the clay gels which have been treated like a homogeneous phase. In fact, as observed by [START_REF] Keller | The pore structure of compacted and partly saturated MX-80 bentonite at different dry densities[END_REF] and [START_REF] Pusch | The microstructure of MX-80 clay with respect to its bulk physical properties under different environmental conditions[END_REF], the clay gels located in the cores of the honeycomb structure present a lower density and a lower mechanical stability than the clay aggregates composing the structure. The rheological behaviour of the clay gels has therefore to be modelled using a variable viscosity (power law, yield stress for fluid flow) [START_REF] Birgersson | Bentonite Erosion. Final Report from Clay Technology[END_REF]. Additionally, a variable density and anisotropic behaviour would enable to characterize more precisely the effective conductivity tensor of clay gels depending on their density.

A perspective of this work, still under development, is to obtain a 3D reconstruction of a hydrated compacted bentonite from a 2D thin section such as the one employed in our numerical computations. This would allow to compute the effective conductivity as well as other transport properties on a microstructure that reflects the morphology and the long-range connectivity of the real porous network [START_REF] Tahmasebi | Reconstruction of 3D porous media using a single thin section[END_REF].

  Figure 1: (a) : Digitalized micrograph of MX-80 clay, taken from (Pusch and Schomburg, 1999). (b) : Finite Element mesh of 1/4 of the microstructure Ω for soft and dense gels, consisting of 271156 elements with 143862 mesh points.

Figure 2 :

 2 Figure 2: Connected mesopores. Fluid velocity field in Case a (soft gels only).

Figure 3 :

 3 Figure 3: Connected mesopores. Fluid velocity field norm and vector in Case b (soft and dense gels).

Figure 4 :

 4 Figure 4: Unconnected mesopores. Distribution of fields Π 1 (a) and Π 2 (b) in the dense gel (porous matrix) domain.

Figure 5 :

 5 Figure 5: No mesopores. Iso-values of Π 1 (a) and Π 2 (b) in the case of clay matrix limited to soft gels.

Figure 6 :

 6 Figure 6: No mesopores. Iso-values of Π 1 (a) and Π 2 (b) in the case of soft and dense gels.

Table 1

 1 Macroscopic hydraulic conductivities computed with unconnected mesopores and no mesopores.

		K transverse	K longitudinal
	Connected mesopores Soft gels	3.319 ×10 -11 -7.735 ×10 -11 ∼ 0
	Connected mesopores Soft + dense gels	1.207 ×10 -9 -2.683 ×10 -9	∼ 0
	Unconnected mesoporess	0.0950 K c	0.1095 K c
	No mesopores Soft gels	0.0215 K c	∼ 0
	No mesopores Soft + dense gels	0.1294 K c	4.9667 ×10 -5 K c

Table 2

 2 Macroscopic hydraulic conductivities computed for Ca-clay (All values are in m/s).

Table 3

 3 Macroscopic hydraulic conductivities computed for N a-clay (All values are in m/s).

		K transverse	K longitudinal
	Unconnected mesopores	1.9 ×10 -11	2.19 ×10 -11
	No mesopores Soft gels	4.3 ×10 -12	∼ 0
	No mesopores Soft + dense gels	2.588 ×10 -11 ∼ 0


Salah Eddine Elayoubi

Younes Khadraoui

Bruno Baynat

Taoufik En-Najjary

Flow level performance evaluation in mobile networks: Analytical modeling and empirical validation

The purpose of this paper is to identify, through empirical validation, the best analytical model for QoS and capacity estimation in mobile data networks. We first present two different sets of models that have been proposed in the literature: the infinite source, Erlang-like models and the finite source, Engset-like models. While the former models are widely developed and adapted to mobile networks, the latter are less used in the mobile context. We thus derive novel analytical models to complete the finite source theory before moving to the comparison of the different models. We make use of network measurements originating from an HSPA network and compare the modeled performance with the observed one, for each of the available analytical models. Our results show that a Processor Sharing analysis that takes as cell capacity the harmonic average of the achievable throughputs and as traffic inputs the volumes generated in each CQI, gives results that are close to the field measurements. As of the finite source models, we observed a good match on some cells and a mismatch on other cells. This is due to the difficulty of extracting traffic parameters from the field. However, finite source models have a better predictive power than infinite source models, as their traffic characteristics can more easily be linked to changes in the network such as the introduction of new services or new devices.

I. INTRODUCTION

With the explosion of data traffic over mobile networks, Quality of Service (QoS) and capacity prediction is becoming a hot topic for operators. Indeed, robust performance evaluation models are needed in order to be able to predict the evolution of QoS and plan capacity upgrade plans. These latter may consist of adding new radio carriers, deploying new sites or even supplementing a site that has 3G with the new 4G technology.

As the objective of mobile operators is to ensure an acceptable user QoS, performance models at flow level are needed. By flow level we mean that the dynamic users' behavior is taken into account. Indeed, in modern data networks, a user initiates a call, transfers a certain amount of data, and leaves the system. The pioneer work of Bonald et al. [START_REF] Bonald | Wireless Downlink Data Channels: User Performance and Cell Dimensioning[END_REF] proposed an analytical model, based on the Processor Sharing (PS) theory, that takes into account the heterogeneity of radio conditions over the cell. In this model, users will accumulate at cell edges as a transfer will take more time if the offered throughput is lower. The authors in [START_REF] Rong | Performance Evaluation of Cellular Networks Offering TV Services[END_REF] extended the model of [START_REF] Bonald | Wireless Downlink Data Channels: User Performance and Cell Dimensioning[END_REF] to the case where multiple classes of service (voice, streaming and data) share the radio resources. Note that these models give a cell capacity that is proportional to the harmonic mean of the throughputs observed over the cell surface, as this harmonic mean gives more weight to positions with lower rates [START_REF] Bonald | How Mobility Impacts the Flow-Level Performance of Wireless Data Systems[END_REF]. Other models like that proposed in [START_REF] Gora | Resource management issues for multi-carrier relay-enhanced systems[END_REF] use an arithmetic mean of throughputs as a measure for the capacity, the main arguments being that user mobility allows observing several radio conditions during a communication, thus preventing user accumulation at cell edge. Recent works have extended Processor sharing models to adaptive streaming traffic [START_REF] Lin | Impact of Chunk Duration on Adaptive Streaming Performance in Mobile Networks[END_REF] and to different scheduling schemes [START_REF] Abbas | Mobility-driven Scheduler for Mobile Networks Carrying Adaptive Streaming Traffic[END_REF], or to take into account mobility of users [START_REF] Baynat | Performance Model for 4G/5G Networks Taking Into Account Intra-and Inter-Cell Mobility of Users[END_REF]. These models have also been extended to the cell coordination schemes in HSPA networks [START_REF] Khlass | Performance Evaluation of Intra-Site Coordination Schemes in HSPA+ Networks[END_REF].

The above mentioned models suppose an infinite source of users: Users arrive from the outside to the network following a Poisson process that is independent from the network state; they are called Erlang-like models, in reference to the original Erlang model proposed by A. K. Erlang [START_REF] Brockmeyer | The Life of A. K. Erlang[END_REF]. Other performance evaluation models for mobile networks are inspired from the Engset model [START_REF] Cohen | The generalized Engset formulae[END_REF] and take the hypothesis of a finite source of users: Users that are physically present in the cell are limited and their activity is described by an ON/OFF process. Examples of such works are [START_REF] Baynat | An efficient analytical model for the dimensioning of wimax networks[END_REF], [START_REF] Baynat | Analytical Models for Dimensioning of OFDMA-based Cellular Networks Carrying VoIP and Best-Effort Traffic[END_REF], that considered that each cell of the network has a finite number of users, each of them generating ON sessions carrying a pre-determined volume of data distributed following an exponential law that is independent from others' activities. Each ON session is then followed by an OFF duration, called reading duration. Note that the available models in the literature suppose that each user may visit all the possible radio conditions during its communication; which is somehow related to a high mobility assumption.

In this paper, we extend these models to the case with low mobility, where users stay with constant radio conditions during a typical communication, thus accumulating at cell edges. Then, in order to identify the best suitable analytical model among the ones of the literature and the new ones developed in the current paper, we perform an empirical validation based on measurements collected from a live HSPA network. The measurement data include radio parameters, indicating the distribution of radio conditions, for each of the considered cells, as well as traffic data. The latter include

A C C E P T E D M

A N U S C R I P T 2 offered traffic volumes per radio condition and typical behavior per user (ON/OFF characteristics). We apply these radio and traffic measurements as inputs to the different analytical methods, and compare the performance results with the Key Performance Indicators (KPI) observed on the field.

The original contributions of this paper are as follows:

• We present and classify the main flow level performance evaluation methods that have been proposed in the literature.

• We extend the Engset-like methods to the case of elastic traffic where users do not change their average radio conditions during a communication.

• We propose a methodology to feed performance evaluation tools with inputs from live networks, in terms of radio conditions and traffic parameters.

• We compare empirically the different analytical models and show the advantages and drawbacks of each of them.

To the best of our knowledge, this is the first work that assesses empirically the performance of flow level models based on field measurements. The remainder of this paper is organized as follows. In Section II, we present the infinite source queuing models. Section III deals with the finite source queuing models and extends the theory to the low mobility case. Section IV shows how the analytical models can be adapted to take as inputs network measurements. Section V compares the outputs of the analytical models to the performance observed in the field and discusses the pros and cons of each method. Finally, Section VI concludes the paper.

II. INFINITE SOURCE QUEUING MODELS

This set of models assumes that there is an infinite source that generates connections in the cell. The theory is well elaborated in the literature, so that we give only an overview of the results that have been obtained in, e.g. [START_REF] Bonald | Wireless Downlink Data Channels: User Performance and Cell Dimensioning[END_REF] and show how that can be applied to an arbitrary number of classes of radio conditions.

A. Radio model

In order to model the capacity of mobile networks, we have to understand how the cell resources are shared among users and how users benefit from obtained resources. We begin by considering round robin scheduling, opportunistic scheduling is considered next. In this case, for a given number of active users, resources (for instance Time Slots in HSPA and Resource Blocks in LTE) are equally divided among users. A user that is alone in the cell will have different bit rates if he is close to the base station, compared to the case where he is far from it, as illustrated in Figure 1. A typical cell of the network can thus be divided into zones of equal radio conditions, or classes, each characterized by an achievable throughput, i.e., a throughput that can be obtained by a user when scheduled by the base station.

Suppose now that the radio conditions are known, obtained analytically, by simulations or from network measurements as will be shown later. We thus assume that the cell can be divided into K zones, each one being associated to a given class i of users, i = 1, ..., K. A user of class i will obtain a throughput C i , if it is alone in the whole cell. These values are not sufficient for estimating the average capacity of the cell. Indeed, the radio interface is a shared broadcast medium, and the traffic dynamics are to be taken into account when considering the capacity, as explained in the following subsection.

B. Processor sharing analysis

We consider here that connection demands arrive to the cell according to a Poisson process of intensity λ connections per second. We assume that a user that carries a new connection demand has a probability p i to be off class i. In addition we suppose that users do not change their radio conditions during their communication and thus remain of the same class during their whole data transfer, i.e., no mobility is considered here. As a result, connection demands of class i arrive to the cell according to a Poisson process of rate λ i = λp i . Each connection (regardless of its class) brings an average amount of data equal to F bits (a file size may follow an arbitrary distribution with average F ), and it ends upon completion of this download. The amount of traffic volume offered to the cell per unit time is thus equal to V = λF (in bit/s). The load generated by users of class i is equal to ρ i = λpiF Ci , and the overall cell load is [START_REF] Bonald | Wireless Downlink Data Channels: User Performance and Cell Dimensioning[END_REF]:

ρ = K i=1 ρ i = λF C (1)
where C corresponds to the harmonic mean of achievable throughput, and can be regarded as the system capacity [START_REF] Bonald | Wireless Downlink Data Channels: User Performance and Cell Dimensioning[END_REF]:

C = 1 K i=1 p i C i (2)
The idea behind this is that, as far as we considered elastic services where the aim is to download a pre-determined amount of data, users at cell edge will have lower data rates and will stay longer in the cell and contribute more to the

PI (n) = ρn (1 -ρ) = V C n (1 - V C ) (3)
the subscript I standing for "infinite source". Knowing the steady-state probabilities, one can derive several performance metrics. For instance, the average number of active flows is calculated by:

n = ∞ n=1 n PI (n) = ρ 1 - ρ (4)
The flow throughput of users of class i is calculated by:

γi = C i (1 -ρ) (5)
and the average flow throughput over the cell is:

γ = C(1 -ρ) (6)

C. Harmonic mean versus arithmetic mean for the cell capacity

Some works in the literature suppose that the cell capacity is given by the arithmetic mean of achievable throughputs Ĉ = K i=1 p i C i , instead of the harmonic mean given in relation 2. Obviously if the probabilities p i involved in averaging Ĉ are the same as those appearing in C, Ĉ will be different from C, and modeling the cell by a PS queue of capacity Ĉ would provide different and thus false results if the aforementioned assumptions (e.g., no mobility of users) are satisfied. However, it is interesting to highlight that an arithmetic mean can provide an equivalent calculation of the cell capacity, if the probabilities involved in the summation have a different interpretation. Let us consider the following expression for Ĉ:

Ĉ = K i=1 p i C i (7)
where p i is not anymore the probability for a user of class i to generate a new connection demand (like for p i ), but the probability that an ongoing transfer in the cell is of class i, i.e., uses throughput C i . p i can alternately be seen as the proportion of ongoing transfers that are of class i. Where p i was only characterized by the ratio λi λ , p i must now be expressed as the ratio ρi ρ , where ρ i = λiF Ci and ρ = K i=1 ρ i . Indeed the more connection demands of class i, the more likely an ongoing transfer is of class i (like in p i ), but in addition, the higher the throughput C i , the faster class i users finish their download and the less likely an ongoing transfer is of class i. Now if we replace in Ĉ the probabilities p i by their expressions, we obtain that Ĉ = C:

Ĉ = K i=1 ρ i C i ρ = K i=1 λ i F ρ = λF K i=1 ρ i = 1 K i=1 p i C i = C (8)
As a conclusion of this subsection, one must be careful when choosing the averaging method, depending on the traffic inputs that are available (volume per radio condition or corresponding load). We will describe the traffic inputs in the numerical evaluation section. And it is of high importance to precisely characterize the probabilities that are involved in the averaging.

D. Introducing opportunistic scheduling gain

The above analysis supposes that the throughput achieved by a user of class i is equal to the throughput when he is alone in the cell divided by the number of active users. This makes the cell capacity C independent of the number of active users. This is a realistic assumption as long as a round robin scheduling is considered. However, opportunistic schedulers like Proportional Fair Scheduler (PFS) are usually activated on the networks and bring a multi-user diversity gain that has to be taken into account in the capacity analysis. The scheduling decision is in this case based on the instantaneous throughput r u (t) of user u at time t (which is zero if the user is not scheduled) and on the average throughput R u (t) of user u at time t, evaluated through the moving average:

R u (t) = (1 -θ)R u (t -1) + θr u (t),
for some fixed time parameter θ (t is discrete and corresponds to scheduling time slots of the user). The PFS schedules the user with the best relative radio conditions, that is with the maximum ratio:

r u (t) R u (t) .
The scheduling gain compared to a blind algorithm like roundrobin depends on the number of users n and on the channel model. If we consider Rayleigh fast fading, the gain G(n) is approximately equal to the harmonic sum 1 + . . . + 1 n for n users [START_REF] Berggren | Asymptotically fair transmission scheduling over fading channels[END_REF]. Note that the gain can be obtained numerically for more realistic channels like PA3 or VA3 channels (as defined in the 3GPP release 12 specifications [START_REF]Study on Further Enhanced Uplink (EUL) enhancements[END_REF]). In this paper, we introduce a scheduling gain calculated following the methodology of [START_REF] Combes | Cross-layer analysis of scheduling gains: Application to LMMSE receivers in frequency-selective Rayleigh-fading channels[END_REF] and that corresponds to a PA3 channel that is usually reported to give realistic results in an urban environment [START_REF] Asplund | How typical is the typical urban channel model? Mobile-based delay spread and orthogonality measurements[END_REF]. The scheduling gains that are used are presented in Table I When opportunistic gain is considered, the system can be modeled as a State Dependent Processor Sharing queue, and the steady state probabilities for the no mobility case will still have a product form given by:

P PFS I (n) = P PFS I (0) ρn n m=1 G(m) (9)
where P PFS I (0) is a normalizing constant. Note that an underlying assumption on users' channels for the product-form expression 9 to hold is that the rate variations are symmetric

A C C E P T E D M A N U S C R I P T

around the mean rate [START_REF] Borst | User-Level Performance of Channel-Aware Scheduling Algorithms in Wireless Data Networks[END_REF]. A similar expression can be found for the high mobility system, by replacing the harmonic average by an arithmetic one in the load expression.

III. FINITE SOURCE QUEUING MODELS

These models suppose that a finite number of users are physically present in the cell. This implies that, unlike the infinite source model, the arrival rate of connections depends on the number of users that are currently being served. Indeed, the more users are in service the less are likely to come. The other key assumption is that, for each user the generated traffic is an ON/OFF process, where the ON period corresponds to an activity session characterized by the average size of the transfered data X on (note that for comparing infinite source and finite source models, we will take X on = F ), while the OFF period corresponds to a reading time characterized by its average duration T off . We consider the same radio model as for the infinite source case, i.e., the cell is characterized by a set of K radio conditions, each one corresponding to a throughput C i .

A. High mobility model

In [START_REF] Baynat | An efficient analytical model for the dimensioning of wimax networks[END_REF] and [START_REF] Baynat | Analytical Models for Dimensioning of OFDMA-based Cellular Networks Carrying VoIP and Best-Effort Traffic[END_REF], authors assume that there is a constant number N of users generating an ON/OFF trafic (as described above), and that, at any time slot, any user that is currently in ON period has a probability p i to be of class i, i.e., to use a throughput C i . In a sense, this corresponds to a high mobility pattern. It has been shown that the cell capacity can be expressed as an arithmetic mean of all achievable throughputs: Ĉ = i p i C i , and the steady state probability of having n ≤ N active users is given by relation [START_REF] Bonald | Network performance analysis[END_REF] [START_REF] Baynat | An efficient analytical model for the dimensioning of wimax networks[END_REF]. Note that, as discussed in Section II-C, the probabilities p i involved in the arithmetic mean of Ĉ does not have the same meaning as those involved in the harmonic mean of C for the infinite source PS queue model (relation ( 2)).

PF (n) = PF (0) X on T off Ĉ n n i=1 (N -i + 1) (10)
the constant PF (0) can be obtained by normalization as:

N n=0 PF (n) = 1
We can finally calculate the cell load as: ρF = 1 -PF (0).

B. No mobility model

We now consider the other extreme case where users cannot change their radio conditions at each time slot, but are rather fixed. To the best of our knowledge, this system has not been modeled in the literature. The system corresponds to the case where there is a fixed number of users for each class i, denoted by N 1 , N 2 , ..., N K . We model this system with a network of K pairs of queues where each pair represents a class of users, be they in an OFF or an ON state, as illustrated in Figure 2.

Proposition 1: The probability of having the distribution n = (n 1 , n 2 , ...n K ) is given by:

P r( n)

= P r[0] n! n 1 !...n K ! K i=1 N i ! (N i -n i )! α ni i ( 11
)
where α i = Xon T off Ci . Proof: The service rate of queue i is equal to:

µ i = C i X on
n i n where n = n 1 + ... + n K is the total number of users. These service rates can be easily shown to verify the symmetric balance equations:

µ i ( n)µ j ( n -e i ) = µ j ( n)µ i ( n -e j ) = n i n j n(n -1)
where e k is a vector with 1 in position k and 0 elsewhere.

The system of 2K queues represented in Figure 2 corresponds thus to a closed Whittle network [START_REF] Bonald | Network performance analysis[END_REF]. It has been shown in [START_REF] Bonald | Network performance analysis[END_REF] that the embedded Markov chain in such a system is reversible and the system has a product-form solution for the steady-state probabilities. These latter can thus be shown, after some calculations, to be equal to that of expression [START_REF] Borst | User-Level Performance of Channel-Aware Scheduling Algorithms in Wireless Data Networks[END_REF].

C. Alternative simple model

The problem of the model with no mobility is that the steady-state probabilities of equation ( 11) cannot be reduced to a single dimension as for equations (3) and [START_REF] Bonald | Network performance analysis[END_REF]. The computation time for such models becomes prohibitive, especially for network dimensioning tools that deal with a large number of cells. In order to cope with this issue, we consider an alternative system that corresponds to a theoretical mobility case, by supposing that users can change their radio condition only after an ON session. This is obviously a simplifying assumption, which makes this model an intermediate case between the high mobility and the no mobility model. We still assume that there is a constant number N of users

A C C E P T E D M A N U S C R I P T 5
generating ON/OFF trafic in the cell, and we define p i as the probability that a user that ends an OFF period is of class i, i.e., uses a throughput C i , during the whole duration of its next ON period. Note that the probabilities p i have, once again, a different interpretation than the probabilities used in the expression of the harmonic mean of C (used in the infinite source PS model), or the arithmetic mean of Ĉ (used in the finite source high mobility model).

This system can be modeled by a network of K + 1 queues, K of them representing the ON states of the K classes, and the remaining one the OFF state, as illustrated in Figure 3. Proposition 2: The probability of having the distribution n = (n 1 , n 2 , ...n K ) is given by:

P r( n) = P r[0] N ! (N -n)! n! n 1 !...n K ! α n1 1 α n2 2 ...α nK K ( 12
)
Proof: The service rate of each queue remains the same as before; the system is thus still a Whittle network [START_REF] Bonald | Network performance analysis[END_REF]. The same analysis as that of Proposition 1 leads to the expression [START_REF] Brockmeyer | The Life of A. K. Erlang[END_REF].

Knowing these steady-state probabilities, we show in the following proposition that the system evolution can be described by the one-dimensional state n = n 1 +...+n K instead of the K-dimensional state n.

Proposition 3: The probability of having n users in the cell can be calculated by:

PF (n) = PF (0) X on T off C n n i=1 (N -i + 1) ( 13
)
where C is the harmonic average of the throughputs over the cell.

Proof: The probability of having n users in the cell is equal to the probability of all possible distributions of n such as n 1 + n 2 + ...n K = n:

P (n) = n|n1+n2+...nK =n P r( n) = P r[0] N ! (N -n)! n|n1+n2+...nK =n n! n 1 !...n K ! α n1 1 α n2 2 ...α nK K
In order to prove equation ( 13), it is sufficient to use the fact that :

n|n1+n2+...nK =n n! n 1 !...n K ! p 1 C 1 n1 ... p K C K nK = K i=1 p i C i n
Note that this expression corresponds to the model with high mobility (equation ( 10)), by replacing the arithmetic average by a harmonic one, and using probabilities p i that have a different interpretation.

The system of Figure 3 can be regarded as an approximation of the system in Figure 2 when the maximal number of users in the two cases verify N = K i=1 N i and N i ≈ p i N . Indeed, as a user during its activity period in Figure 3 does not change its radio class, the phenomenon of accumulation of users with bad radio conditions, characterizing the system of Figure 2 will still happen. In order to illustrate this, we plot in Figure 4 the cell load obtained for different values of X on for the approximate and exact finite source models with static users (equations ( 13) and ( 11)), taking the same simple example of the previous section, with T off = 40 seconds and N = 10 users in the cell. We observe that the harmonic average method is a good approximation for the system with no mobility. We thus use, in the remainder of this paper, the harmonic average method for modeling the no mobility case.

D. Impact of opportunistic scheduling

As explained for the infinite source models, multi-user diversity gains are to be taken into account in the performance calculations.

Proposition 4: The steady-state expressions ( 10) and ( 13) can be extended for considering a scheduling gain G(n) that depends on the number of active users, as follows:

P PFS F (n) = P PFS F (0) λ u X on Ĉ n 1 n m=1 G(m) n i=1 (N -i+1) (14)

A C C E P T E D M

A N U S C R I P T 6 and

P PFS F (n) = P PFS F (0) λ u X on C n 1 n m=1 G(m) n i=1 (N -i+1) (15
) Proof: The Whittle network condition still holds when the service rate of users of class i is equal to:

µ i = C i X on n i n G(n)
Expressions ( 10) and ( 13) can be derived as for the case of a constant capacity.

IV. EMPIRICAL DATA DESCRIPTION

The multiplicity of analytical models that can be used for performance evaluation in cellular networks shows the need for an empirical validation and comparison of these models. In this section, we describe the measurements that have been collected from the field for being used in the validation process. We consider a live HSPA network in a major European city, focus on the downlink and make use of two sources of data:

• Operation and Maintenance Center (OMC) counters:

these counters aggregate the information reported by user equipments to base stations, in addition to the data measured directly by the base stations. They include traffic volume information, radio conditions and performance metrics (loads, number of active users, etc).

• Probes data: Detailed traffic information that describe the behavior of users (number of connections per user, traffic volume per connection, etc.) cannot be obtained from OMC counters, but rather from traffic probes installed on the interface between the access and the core network.

Our dataset consists of two captures of data call records of one hour each, the first one between 9 am and 10 am and the second trace collected between 7 pm and 8 pm. For each radio connection, the probe collects several pieces of information, namely: customer id (assigned anonymously by the probe, the probe keeps the same id inside the area covered by the RNC), the start and end time of the connection, DL/UL volumes in bytes, and cell identifier (assigned anonymously by the probe).

A. Radio measurements

We first begin by describing the radio-related parameters used in the validation process.

• CQI distribution: The main signal quality measure in HSPA is the so-called Channel Quality Indicator (CQI). This measurement is a quantified Signal -to-Interferenceplus-Noise-Ratio (SINR) reported from UEs to their serving base station, available later at the OMC level. Using these field measurements, a CQI distribution can be constructed. An example CQI distribution obtained from the field is given in Figure 5.

• HSPA Device mix: This distribution indicates the proportion of HSPA device categories in terms of traffic volume. The construction of this distribution is similar to that stated for CQI distribution, with the appropriate field counters. Note that HSPA devices that are currently present on the networks range from category 6 (16 QAM modulation with a peak throughput of 3.6 Mbit/s) to Category 24 (dual cell with 64 QAM and a peak throughput of 42.2 MBit/s).

B. Traffic measurements

The simplest traffic measure to obtain is the overall traffic intensity V in bit/s. It can be obtained by dividing the traffic volume, a KPI provided by OMC counters, by the observation time.

The other traffic measures that are used in this paper, especially for the finite source models, cannot be obtained from OMC counters, but from probes. Before estimating the average session volume X on and average time between sessions T off we need to take into account some characteristics of the new generation of smartphones. Indeed, new smartphones are designed to be always connected to Internet in order to receive notifications for the applications installed on them like social networks, email, news, etc. To receive notifications, smartphones either keep a continuous TCP connection with the servers by regularly sending keep-alive packets, or regularly initiate a new TCP connection toward the servers to check for new notifications. These background activities generate a huge number of very short radio connections with negligible amount of data traffic. Separating background activity from regular activity is not an easy problem. In this work, we assume all radio connections with less the 2 Kbytes (less than two MTU TCP packets) as background activity. These connections represent 62% of total number of radio connections but only 0.001% of transmitted volume. Keeping these connections while doing the estimation of X on and T off lead to a ridiculously small values that are not representative of the actual connections generating the observed traffic volumes. Therefore, we drop all these connections from our dataset.

In order to make the calibration process easier, we suppose that X on and T off are characteristics of customers behavior and are not dependent on the network cells. This of course a simplification (that will be relaxed in Section V-B), as user behavior is known to actually depend on the cell [START_REF] Trestian | Measuring Serendipity: Connecting People, Locations and Interests in a Mobile 3G Network[END_REF]. Therefore, we first estimate these values over the whole trace after cleaning out the background activity. We thus obtain, for each operating hour, specific values of X on and T off . We will check the validity of this assumption afterwards. We perform the following estimations:

• In order to determine if two consecutive packets having the same destination belong to the same downloading ON period, or to two consecutive ON periods separated by a reading OFF period, we considered a constant threshold (empirically set to 2s) under which the two packets are assumed to belong to the same ON period, and above which they are assumed to belong to two consecutive ON periods.

• X on is estimated as the average of the volumes of all the radio connections in our dataset.

• For each customer in our dataset (we have more than 80 thousands customers per trace) we estimate the median interconnection time (MIT), and T off is estimated as the median of all MITs of all customers.

• For each cell, the number of users physically present in the cell N is estimated as the total number of customers connected to this cell during our observation. We make use of two estimates: the first takes into account all the unique users observed during the considered hour on the cell; this gives an absolute maximum of the number of physically present users. The other estimate considers only users that generated at least two connections during the considered hour, in order to be sure that these users stayed in the cell for a significant amount of time.

C. Performance metrics

In addition to the traffic inputs and radio parameters, OMC data give some performance metrics, for instance:

• Scheduler load: this OMC metric corresponds to the proportion of time slots that are used for data transmission.

• Average number of active users: this OMC metric gives the average number of users for which the base station has some data to transmit (that are scheduled or wait to be scheduled).

• Average user throughput: this is not a metric that is directly given by counters, but can be estimated as the ratio between the traffic volume and the average number of active users.

D. Taking into account voice traffic

The models presented in this paper are specific to elastic data traffic. However, HSPA networks usually share the same spectrum resources with circuit switched services (R99 voice service) and some cells have a significant volume of voice traffic. As the circuit-switched traffic has priority in power resource allocation, the remaining power resource is assigned to HSPA, which reduces the throughput of data users compared to a carrier that is dedicated to HSPA. In order to estimate the amount of power consumed by voice users, we make use of multi-Erlang theory, as in [START_REF] Baroudy | HSUPA/HSDPA systems: capacity and dimensioning[END_REF]. Indeed, if we keep the same classification of the cell into zones of equal radio conditions, let P v i be the power consumed by a voice user at position i, the amount of power consumed by voice users when there is a vector of n v = (n v 1 , ..., n v K ) voice users in the cell is equal to:

P v ( n v ) = K i=1 n v i P v i ≤ P max -P cont
where P max is the maximal transmitting power of the base station and P cont is the power reserved for control channels. Multi-Erlang analysis (using Kaufman-Roberts algorithm [17][20]) allows computing the blocking rate for voice users and their average consumed power Pv . The remaining available power for HSPA traffic is thus P max -P cont -Pv . This has an impact on the throughput that can be achieved by HSPA users (the achievable throughput is lower than the case of a dedicated carrier for HSPA), and this reduction can be taken into account by a shift of the CQI distribution by a number of dBs equivalent to the power reduction.

In order to take into account the impact of voice traffic on the performance, we make use of the following OMC counters:

• The voice traffic volume (in Erlang).

• The maximal cell power.

• The common channels power. As of the powers consumed per voice call in different positions in the cells, it is not given by any field counter. We thus make use of the simulation results reported in [START_REF] Baroudy | HSUPA/HSDPA systems: capacity and dimensioning[END_REF] for estimating them.

E. Adaptation of measurements to analytical models

After describing the data that is available in the field, we now describe how these measurements can be used as inputs for the analytical models of Sections II and III.

First, we note that the analytical models described in Section II take as input a set of achievable throughputs, {C i }, with their associated weights p i . However the inputs collected from measurements give a distribution of CQIs, as explained above. In order to obtain the achievable throughputs, we make use of Transport Block Size (TBS) tables provided by the 3GPP, that associate each CQI to a volume of traffic per TTI (or equivalently a throughput in bit/s), for each device category [24]. Combining the device mix observed on the field with the CQI mix, the achievable throughput distribution can be obtained. As of the p i 's, they correspond to the weight associated to each CQI and represent the percentage of users that see the different CQIs. They correspond thus to a proportion of offered volume per radio condition.

Table II shows the measurement inputs used by the different analytical models and the outputs that can be compared to the KPIs observed on the field. Both infinite and finite source models take as input the radio parameters (CQI and device mix, powers). However, only one HSPA traffic parameter -the traffic volume -is needed for the infinite source models, while three traffic parameters are needed for the finite source models (average session volume, average time between sessions and number of users that are present in the cell). Some other KPIs observed on the field (average cell load, number of active users, average user throughput) are used to compare the models outputs to the field.

V. VALIDATION RESULTS AND DISCUSSIONS

A. Infinite source models

We begin by presenting the results obtained with infinite source models. The models are tested on a set of 7 cells, selected so that they represent the typical traffic load values observed in the network during the measurement hours. Figures 6 and 7 present the scheduler loads and the numbers of active users, respectively. The x axis of these figures represents a fictitious cell ID, where cells are sorted following their increasing observed loads. Three values are plotted for each of the cells corresponding to the field measurements, the result obtained using a harmonic average of throughputs as the cell capacity and the result obtained using an arithmetic average (with the same weighting probabilities, as suggested in [START_REF] Gora | Resource management issues for multi-carrier relay-enhanced systems[END_REF]). The first observation is that the harmonic average gives in general acceptable results that are close to the field measurements, for both loads and average number of active users. The second observation is that the arithmetic average method gives optimistic results (low loads and small numbers of active users). We explain this by the fact that the majority of the traffic in the considered HSPA network (an in mobile networks However, there is still some gap between the field observations and the analytical results with a harmonic average; we can give some insight into the possible causes of this inaccuracy:

• Mobility: Even if the majority of traffic is generated by non moving users, some of the traffic is generated by users in mobility. Neither the harmonic average nor the arithmetic one are able to describe the performance. However, the results are closer to the harmonic average method as there are more static users than mobile ones.

• Inaccuracies in the inputs: Even if we tried to stick the best to the field measurements for the inputs of the analytical models, there are still some inputs values that we obtained from simulations as they cannot be directly observed on the field. This includes scheduling gains (Table I) and the power consumed by a voice user.

• Equipment specificities: Our models consider generic radio resource management (scheduling, throughput association for each of the reported CQIs, etc.) that does not take into account the specificities of the different vendor implementations.

• Devices with advanced receivers and receive diversity:

The device mix reported by the OMC counters does not include the information about receivers and receive diversity. We thus make a conservative assumption of standard receivers with no receive diversity, which tends to give pessimistic results.

B. Finite source models

We now move to the finite source model. We consider the model presented in equation [START_REF] Cohen | The generalized Engset formulae[END_REF] for two values of the number of users N : the upper bound (all users) and the number of users that made at least two data connections during the observation period. Figures 8 and9 present, respectively, the scheduler loads and the numbers of active users, and follow the same presentation as before.

We first observe that the gap between results obtained when considering all users versus the performance when filtering on users with multiple connections is high: Taking all users, even those that make only one connection during the observation hour, increases the calculated cell load. The second observation is that when there is a match between the field observation and the prediction, the filtered case is the closest one; this is the case of cells 1, 4 and 6. For the other cells, we observe that both values (all users and filtered connections) give optimistic results compared to the field (lower loads and lower number of users). This may result from a sub-estimation of the other traffic parameters: the volume per connection X on and the OFF duration T off . We indeed based the X on and T off calculations on the RNC basis and not on a cell per cell basis, arguing that the behavior of users is not cell-dependent and that the differentiating parameter between cells is the number of users.

In order to investigate this issue of user's behavior that is different in some cells, we focus on cell 7 and extract the traffic characteristics of users that are connected to this cell (the results are not shown in Figures). We obtain an average connection size X on of 315 Kbytes, compared to the RNC average of 243 Kbytes per connection, while the T off value remains almost stable. Applying this new traffic characteristic to this cell, we obtain a load of 50%, very close to the 56% observed on the field (this is equivalent to the result obtained with the infinite source model using the harmonic average). This experiment illustrates the fact that the finite source models, if tuned on a cell per cell basis, can give good results that fit well to the field.

These results illustrate the difficulty of calibrating finite source models with field measurements. Indeed, while only the traffic volume V is needed to be extracted on a cell-by-cell basis in infinite source models, three parameters are needed in finite source models: X on , T off and N . Note however that analytical models with finite sources are only sensitive to the ratio Xon T off .

VI. CONCLUSION

The objective of this paper is twofold: modeling the cell capacity and comparing the two main sets of models for performance evaluation in mobile networks. We followed an empirical approach that consists of feeding the analytical As of the first objective, we considered the two definitions of capacity in the literature that take into account the heterogeneity of radio conditions over the cell. We namely consider as capacity the arithmetic average of achievable throughputs, versus their harmonic average, and highlight the fact that, before comparing them, one must precisely characterize the probabilities that are involved in the averaging. We show that the harmonic average corresponds to the static case, while the arithmetic average is more related to the high mobility case. Our experiments show that the analytical results match better to field observations when the capacity is considered as the harmonic average of the achievable throughputs, as recommended in [START_REF] Bonald | Wireless Downlink Data Channels: User Performance and Cell Dimensioning[END_REF]. This is due to the fact that a large proportion of the traffic in mobile networks is generated by indoor or fixed users. For example, [START_REF]G PPP association, 5G-PPP use cases and performance evaluation models[END_REF], Table 13, states that 80% of traffic in dense urban environments is generated by indoor users.

For the model itself, we studied the infinite source models where the arrival rate of connections is supposed to be independent of the network state, and the finite source models where the number of users that are physically present in the cell is limited so that the arrival rate reduces when the number of active users increases. We presented the models that are available in the literature for infinite source models and derived novel finite source queuing models that correspond to the case of a low mobility. We showed that finite source models are pretty difficult to calibrate with field measurements, due to a larger number of needed traffic inputs on a cell per cell basis, while infinite source models are simpler to calibrate because they only take into consideration the traffic volume. However, finite source models are able to capture the impact of user behavior and are useful, for instance, in capacity/technology planning based on marketing predictions. In this context, finite source models are very useful for predicting the impact on the performance of the introduction of new services/devices, as these latter usually introduce changes in user behaviors (activity patterns, connection volumes, etc.).

Fig. 1 .

 1 Fig. 1. a) Example of simulated achievable throughput over an HSPA cell: only one user is in the cell but is simulated at different positions in the cell. b) Cell decomposition into zones of equal radio conditions.

  Based on this load expression, one can model the cell as a Processor Sharing (PS) queue of capacity C, and the steady-state probability of having n active users in the cell are calculated by:

Fig. 2 .

 2 Fig. 2. No mobility system: each tandem of queues represents users of a given class. When a user achieves his service (the ON period) he returns to the OFF mode until the next activity period.

Fig. 3 .

 3 Fig.3. Athernative system with a theoretical low mobility pattern: When a user achieves his OFF period, he choses following probability p i to join queue i during his ON session; he then returns to the OFF queue until his next activity period.

Fig. 4 .

 4 Fig. 4. Cell load vs. Xon

Fig. 5 .

 5 Fig. 5. An example field CQI distribution

  . avg. PS: Arith. avg.

Fig. 6 .

 6 Fig. 6. Validation results for infinite source models: cell loads

Fig. 7 .

 7 Fig. 7. Validation results for infinite source models: number of active users

Fig. 8 .

 8 Fig. 8. Validation results for finite source models: cell loads (the same cells as before are used).

Fig. 9 .

 9 Fig. 9. Validation results for finite source models: number of active users

* This work has been carried out in the framework of IDEFIX project, funded by the ANR under the contract number ANR-13-INFR-0006.