Aliaksei Pazdniakou
email: aliaksei.pazdniakou@univ-lorraine.fr

Anne-Julie Tinet
email: anne-julie.tinet@univ-lorraine.fr

Fabrice Golfier
email: fabrice.golfier@ensg.univ-lorraine.fr

Kassem Kalo

Stephane Gaboreau

Patrick Gaire

Numerical efficiency assessment of the lattice Boltzmann model for digital nano-porous rock applications

Keywords: nano-porous media, Lattice Boltzmann Method (LBM), Klinkenberg effect

Study of transport properties of low permeable sedimentary clay formations is of growing interest in exploration of non-conventional resources, and for long term radioactive waste repositories where these formations play the role of a natural shield which prevents the spreading of radioactive materials in the environment. Recent advances in imaging techniques allow to obtain 3D volume of compacted clay sample with pixel size of order of few nano-meters which enables us to explore transport properties of these media by means of numerical simulation. In our study, the Lattice Boltzmann Model (LBM) is applied to simulate fluid flow inside the 3D image of a compacted illite sample to explore its transport properties. The impact of different image segmentation methods applied to treat the raw data on the numerical results is addressed. Because of the nanometric pore size, various physical phenomena may influence fluid flow. In this study, we consider the influence of gas slippage effect (the Klinkenberg effect) on the permeability. The numerically obtained results are compared with experimental data obtained for a reconstructed compacted illite sample. If the original geometry is not strongly modified by the applied image segmentation method, the numerically calculated Klinkenberg factor follows well existing trends and provides predictions close to the experimental values.

Introduction

The main objective of this study is to improve the understanding of one-phase flow in low permeable porous media by means of numerical simulations performed at the pore scale on 3D images of clay samples. This objective is closely related with the increasing interest in low permeable sedimentary formations which are projected to be used as host site for long term radioactive waste repository [START_REF] Andra Dossier | Evaluation of the feasibility of a geological repository in an argillaceous formation[END_REF]. Characterized by low permeability, these media can serve as a natural barrier against radioactive materials. During the storage phase, indeed, a certain amount of gas (mostly hydrogen) will be released due to corrosion of iron components of the facility. Increased gas pressure can cause gas migration through the highly water saturated sedimentary material and even affect the geometry of the pore space modifying transport properties of the clay. Therefore, a good understanding of this complex process is of high interest.

Depending on hydrodynamic conditions, various transport mechanisms participating to gas migration can be in action, separately or together [START_REF] Marschall | Characterisation of Gas Transport Properties of the Opalinus Clay, a Potential Host Rock Formation for Radioactive Waste Disposal[END_REF][START_REF] Cuss | Experimental observations of mechanical dilatation at the onset of gas flow in Callovo-Oxfordian claystone[END_REF]. Thus, gas diffusion can be favoured, and the dissolved gas will be transported by advection and diffusion [START_REF] Bardelli | Hydrogen uptake and diffusion in Callovo-Oxfordian clay rock nuclear waste disposal technology[END_REF]. Conversely, for sufficiently high pressures, gas percolates through the water saturated clay and is transported by continuous visco-capillary flow or small bubble migration. The excess of gas pressure can lead to the dilatation of existing and creation of new pathways [START_REF] Cuss | Experimental observations of mechanical dilatation at the onset of gas flow in Callovo-Oxfordian claystone[END_REF][START_REF] Harrington | Gas flow in Callovo-Oxfordian claystone (COx): results from laboratory and field-scale measurements[END_REF][START_REF] Davy | Migration de gaz au travers de l'argilite de Bure: mise en evidence d'un passage discontinu[END_REF][START_REF] Pazdniakou | Migration of gas in water-saturated clays by coupled hydraulicmechanical model[END_REF]. Due to nano-metric pore size, however, specific hydrodynamic features need to be considered.

The clayey material has a complex porous structure with the total porosity ranging from 10% to 40% (up to 50% for bentonite), most of which is associated with pores smaller than 100 nm [START_REF] Gaboreau | Optimization of pore-network characterization of a compacted clay material by TEM and FIB/SEM imaging[END_REF][START_REF] Marschall | Characterisation of Gas Transport Properties of the Opalinus Clay, a Potential Host Rock Formation for Radioactive Waste Disposal[END_REF]. If gas pressure is not sufficiently high, the mean free path of gas molecules λ can be significant when compared to the pore size H causing the gas slippage effect. The intensity of gas rarefaction effect is characterized by the Knudsen number Kn Kn = λ H

The continuous flow regime corresponds to Kn < 0.001; the slip flow is associated with 0.001 < Kn < 0.1; the transitional flow regime corresponds to 0.1 < Kn < 10, and the free molecular flow is related to Kn > 10 [START_REF] Mitra | Microfluidics and Nanofluidics Handbook, Two Volume Set[END_REF][START_REF] Guo | Lattice Boltzmann Method and its Applications in Engineering[END_REF]. For the description of slip flow, the Navier-Stokes equations still can be used together with a slip boundary condition. At the macroscopic scale, the gas slippage effect is manifested in the increased gas permeability K, which is known as the Klinkenberg effect [START_REF] Klinkenberg | The permeability of porous media to liquid and gases, American Petroleum Institute, Drilling and Productions practices[END_REF] and historically expressed as

K = K 0 1 + b P m (2)
where K 0 is the absolute permeability, b is the Klinkenberg factor related to the medium properties, and P m is the mean gas pressure.

Recent advances [START_REF] Gaboreau | Optimization of pore-network characterization of a compacted clay material by TEM and FIB/SEM imaging[END_REF][START_REF] Song | Multi-scale pore structure of COx claystone: Towards the prediction of fluid transport[END_REF] in digital rock imaging improved the understanding of the complex hierarchical pore space structure of sedimentary clays. In absence of micrometric fractures, most of the interconnected pores are of nanometric size. The access to the high resolution digital images of real clay samples obtained by FIB-SEM (Focused Ion Beam Scanning Electron Microscope) [START_REF] Gaboreau | Optimization of pore-network characterization of a compacted clay material by TEM and FIB/SEM imaging[END_REF], or at the lower scale by Transmission Electron Microscopy (TEM) if the meso-pores are not connected [START_REF] Davy | Three-scale analysis of the permeability of a natural shale[END_REF], provides the opportunity to perform numerical simulations of fluid flows, in order to better understand the transport properties of these low permeable porous media. Various numerical methods can be employed to calculate capillary pressure, absolute and relative permeabilities. In our study, we focus on using the lattice Boltzmann method [START_REF] Chai | Gas Flow Through Square Arrays of Circular Cylinders with Klinkenberg Effect: A Lattice Boltzmann Study[END_REF]. Originally derived as a generalization of lattice gas automata (LGA), the LBM can be also obtained by discretizing the continuous Boltzmann equation in space, time, and particle velocity space. Unlike the classical Computational Fluid Dynamics (CFD) methods, the LBM describes fluid state by the probabilistic particle distribution function. The model was successfully applied to simulate single-and multiphase flow inside the 3D images of porous media [START_REF] Pazdniakou | Dynamic permeability of porous media by the lattice Boltzmann method[END_REF][START_REF] Dymitrowska | Two-phase flow pore-size simulations in Opalinus clay by the lattice Boltzmann method[END_REF][START_REF] Benioug | An immersed-boundary lattice Boltzmann model of biofilm growth in porous media[END_REF][START_REF] Ben Hadj Hassine | Gas Migration in Highly Water-Saturated Opalinus Clay Microfractures Using a Two-Phase TRT LBM[END_REF]. The gas slippage effect was introduced later by adopting an appropriate solid-fluid boundary condition. Different formulations have been derived for this purpose based on the pure diffuse reflection [START_REF] Ansumali | Kinetic boundary conditions in the lattice Boltzmann method[END_REF], its combination with the specular reflection [START_REF] Tang | Lattice Boltzmann method for gaseous microflows using kinetic theory boundary conditions[END_REF] or with the bounce-back boundary condition [START_REF] Chai | Gas Flow Through Square Arrays of Circular Cylinders with Klinkenberg Effect: A Lattice Boltzmann Study[END_REF], and a combination of the bounce-back with the specular reflection [START_REF] Succi | Mesoscopic modeling of slop motion at fluid-solid interfaces with heterogeneous catalysis[END_REF].

Up to now, however, the use of these boundary conditions is limited to simplified and well discretized geometries (e.g., capillary tubes) with the notable exception of the recent application by [START_REF] Zhao | Study of Gas Flow Characteristics in Tight Porous Media with a Microscale Lattice Boltzmann Model[END_REF] to 3D digital rocks but without comparison to experimental permeability data. However, challenges still remain to gain reliability in predicting transport properties of such nanoporous materials using Digital Rock Physics (DRP) measurement. First, in spite of the increasing technological advances, it is still unclear whether the spatial resolution available by 3D microscopic imaging techniques (X-ray tomography, FIB-SEM, BIB-SEM (Broad Ion Beam Scanning Electron Microscope),...) is sufficient to properly capture rock surface features (pore curvature, roughness) that may impact the slippage effect for computer-based simulations. In other words, concerns may be raised whether the filtering process conducted for measuring integrated properties such as permeability possesses sufficient imaging-based microscale information in the presence of slip flow. In adddition, even if the solution of LBM model with slip flow boundary condition has been validated rigorously against fluid flow solution through a parallel-wall uniform fracture or Poiseuille tube flow, applicability to 3D complex pore geometry with non-uniforme pore throats (i.e, with spatially variable Knudsen number) needs to be assessed. The present work aims at addressing these issues.

Our ultimate objective is to explore the capability of the LBM to predict the transport properties of real nano-porous media and confront them to experimental data. An artificial clay material, initially imaged in 3D by [START_REF] Gaboreau | Optimization of pore-network characterization of a compacted clay material by TEM and FIB/SEM imaging[END_REF] will be considered for this purpose. It must be noted that, with the present computational resources, it is not possible to treat the whole image of the sample accounting for all pore scales from micro-to nanometric. Therefore, we are going to work with a selected sub-sample insuring the pore space connectivity. The sensitivity of the numerical methods to the geometry resolution is also investigated. To the best of our knowledge, this is the very first attempt to tackle this complex issue.

Numerical model

The Navier-Stokes equations describing the flow inside the pore-space can be numerically solved by various computational fluid dynamics (CFD) methods such as finite volume method (FVM) or finite element method (FEM). Recently, non-conventional numerical methods such as LBM received a lot of attention for single-and multi-phase flow modelling in porous media [START_REF] Pazdniakou | Dynamic permeability of porous media by the lattice Boltzmann method[END_REF][START_REF] Dymitrowska | Two-phase flow pore-size simulations in Opalinus clay by the lattice Boltzmann method[END_REF]. Indeed, compared to conventional numerical methods, they allow to cope naturally with fluid-fluid interface and complex fluid-solid boundary conditions. Also, various physical phenomena can be relatively easily incorporated in the basic LBM. Finally, possessing a high level of intrinsic parallelism, the LBM can easily profit from parallel computing on CPU and GPU clusters, which is important when dealing with prohibitively large samples [START_REF] Alpak | A distributed parallel multiplerelaxation-time lattice Boltzmann method on general-purpose graphics processing units for the rapid and scalable computation of absolute permeability from high-resolution 3D micro-CT images[END_REF].

2.1. Lattice Boltzmann Method.

The D3Q19 TRT LBM

The LBM originally appeared as a generalization of LGA in an attempt to overcome some of their difficulties. It was first introduced for a hexagonal lattice with six discrete velocities (D2Q6) in 2D space and a face-centered hypercube (FHCH) lattice in 4D space [START_REF] Mcnamara | Use of the Boltzmann equation to simulate lattice-gas automata[END_REF]. Later, it was demonstrated that LBM can be derived from the continuous Boltzmann equation by applying a proper discretization [START_REF] He | A priori derivation of the lattice Boltzmann equation[END_REF][START_REF] He | Theory of the lattice Boltzmann method: From the Boltzmann equation to the lattice Boltzmann equation[END_REF]. The method rapidly extended to various lattices, collision operators, and boundary conditions. Now, it is applied to simulate a wide variety of fluid flow phenomena.

In our study, a regular cubic lattice in 3D space with 19 discrete velocities (D3Q19) is adopted. In LBM, the state of the fluid at each lattice point r at a given time step t is defined by the particle distribution function f = (f i (r, t), i = 1, ..., Q), whose components are associated with discrete velocities c i . For a given time moment t each f i (r, t) defines the probability to find a particle with the discrete velocity c i at the location r. The 19 discrete velocities for D3Q19 are given by

c i =      (0, 0, 0), i = 1 (±c, 0, 0), (0, ±c, 0), (0, 0, ±c) i = 2, ..., 7 (±c, ±c, 0), (±c, 0, ±c), (0, ±c, ±c) i = 8, ..., 19 (3)
where the velocity c = δx δt is defined as a ratio of the lattice step δx and the time step δt. Both the lattice and the time steps are set to 1 for simplicity. The evolution of the particle distribution function is described by

f i (r + c i δt, t + δt) = f i (r, t) -Ω i (f, f eq ) + F i δt (4)
where Ω i (f, f eq ) is the collision operator, F i is the discretized body force, and f eq is the equilibrium distribution. The algorithm can be easily split into two steps, i.e., the calculation of the post-collision distributions in the right part of (4), and the propagation of the post-collision distributions according to the discrete velocities directions. The collision operator Ω i (f, f eq ) describes the evolution of the particle distribution function due to collision between the particles. The most simple collision operator is the Bhatnagar-Gross-Krook (BGK) collision operator [START_REF] Bhatnagar | A model for collision processes in Gases. I. Small amplitude processes in charged and neutral one-component systems[END_REF] which relaxes the particle distribution function to its equilibrium with a single relaxation time. When used with the bounce-back boundary condition, the well known drawback of this operator is the viscosity dependent location of the boundary. The solution to this problem is the multirelaxation time collision operator (MRT) which relaxes several fluid moments to their equilibrium states with individual relaxation times [START_REF] Humières | Multiple-relaxation-time lattice Boltzmann models in three dimensions[END_REF]

Ω i (f, f eq ) = M -1 • S • M(f -f eq ) ( 5
)
where M is the transformation matrix, and S is the diagonal relaxation matrix. The transformation matrix M projects the particle distribution function onto the moment space giving the following moments [START_REF] Humières | Multiple-relaxation-time lattice Boltzmann models in three dimensions[END_REF] M • f = (ρ, e, ε, J x , q x , J y , q y , J z , q z , 3p xx , 3π xx , p ww , π ww ,

p xy , p yz , p xz , m x , m y , m z ) ( 6
)
where ρ is the fluid density, and J x , J y , J z are the fluid momentum components; e is related to the kinetic energy and ε to the kinetic energy square; q x , q y , q z correspond to energy flux, and 3p xx , p ww , p xy , p yz , p xz to the components of the viscous stress tensor; there are also two fourth order moments 3π xx , π ww and three third order moments m x , m y , m z . For the D3Q19 model, the M matrix can be found in [START_REF] Humières | Multiple-relaxation-time lattice Boltzmann models in three dimensions[END_REF]. The diagonal matrix S defines the relaxation rates for different fluid moments [START_REF] Humières | Multiple-relaxation-time lattice Boltzmann models in three dimensions[END_REF] S = diag(0, s e , s ε , 0, s q , 0, s q , 0, s q , s ν , s π , s ν , s π , s ν , s ν , s ν ,

s m , s m , s m ) (7)
The relaxation rates are positive and must be smaller than 2 for numerical stability.

The dimensionless relaxation times are the inverse of the relaxation rates, for example τ ν = 1 sν . When using the MRT collision operator, the kinematic ν and bulk ζ fluid viscosity can be prescribed separately [START_REF] Humières | Multiple-relaxation-time lattice Boltzmann models in three dimensions[END_REF]

ν = c 2 s ( 1 s ν - 1 2 )δ t (8a) ζ = ( 5 9 -c 2 s )( 1 s e - 1 2 )δ t (8b)
where c s = c √ 3 is the speed of sound. Relaxation times can be tuned in order to improve the precision and numerical stability [START_REF] Humières | Multiple-relaxation-time lattice Boltzmann models in three dimensions[END_REF]. The most popular choice is the two relaxation times (TRT) model with two different relaxation times τ ν and τ q , which relaxes even and odd order fluid moments with different rates [START_REF] Ginzburg | Multireflection boundary conditions for lattice Boltzmann models[END_REF]

s e = s ε = s π = s ν , s m = s q = 8 2 -s ν 8 -s ν (9)
The collision operator relaxes the particle distribution function towards its equilibrium state. In LBM, the low Mach number expansion of the Maxwellian distribution is usually used

f eq i (ρ, u) = w i ρ[1 + 1 c 2 s (c i • u) + 1 2c 4 s (c i • u) 2 - 1 2c 2 s |u| 2 ] ( 10
)
where ρ is the fluid density, u is the fluid velocity, and w i are the lattice weights. This expansion permits to recover the Navier-Stokes equation at the macroscopic level by using the Chapman-Enskog decomposition. The lattice weights w i in [START_REF] Chai | Lattice Boltzmann simulation of surface roughness effect on gaseous flow in a microchannel[END_REF] are given by

w i =    1 3 , c i = 0; 1 18 , c i = 1; 1 36 , c i = √ 2.
(11)

Macroscopic variables

The fluid density ρ and the fluid velocity u can be obtained as moments of the distribution function

ρ = Q i=1 f i = Q i=1 f eq i (12a) ρu = Q i=1 f c i = Q i=1 f eq i c i (12b)
In presence of the body force ρa, the fluid momentum is redefined as [START_REF] Ladd | Lattice-Boltzmann simulations of particle-fluid suspensions[END_REF] ρu

= Q i=1 f c i + ρa 2 δt ( 13
)
in calculation of the nonlinear terms of the equilibrium distribution function as well as the output fluid velocity in order to improve the precision. The pressure is provided by the equation of state of the ideal gas

p = c 2 s ρ (14)

Forcing

The body force ρa must be discretized as F = (F i (r, t), i = 1, ..., Q) to be incorporated in LBM. The most common way is to use the following discretization [START_REF] Luo | Theory of the lattice Boltzmann method: Lattice Boltzmann models for nonideal gases[END_REF]

F i = w i ρ c i -u c 2 s + c i • u c 4 s • a (15)
This discretization allows to use the body force which is not constant in space preserving the Galilean invariance [START_REF] Luo | Theory of the lattice Boltzmann method: Lattice Boltzmann models for nonideal gases[END_REF]. In this paper following [START_REF] Guo | Lattice Boltzmann Method and its Applications in Engineering[END_REF], the discretized body force F given by ( 15) is further relaxed to F as

F = M -1 • I - 1 2 S • M • F (16)
before the incorporation.

No slip boundary condition

The interaction between the fluid and solid phase is described by the boundary conditions. The boundary condition must supply the unknown incoming distribution functions for the fluid nodes neighbouring the solid nodes. The most simple and widespread boundary condition is the bounce-back [START_REF] Ginzbourg | Boundary flow condition analysis for the three-dimensional lattice Boltzmann model[END_REF]

f i (r, t + δt) = fi (r, t) (17)
where fi (r, t) is the post-collision distribution corresponding to the discrete velocity c i = -c i . According to the bounce back rule, the post-collision distribution is reflected from the solid node in the opposite direction. Combined with the TRT collision operator, this rule provides the second order no-slip boundary condition [START_REF] Ginzburg | Multireflection boundary conditions for lattice Boltzmann models[END_REF]. The no-slip surface is located in the middle between the fluid and solid lattice nodes. It is possible to change this position by using various interpolation based boundary conditions [START_REF] Chun | Interpolated boundary condition for lattice Boltzmann simulations of flows in narrow gaps[END_REF][START_REF] Ginzburg | Multireflection boundary conditions for lattice Boltzmann models[END_REF][START_REF] Pan | An evaluation of lattice Boltzmann schemes for porous medium flow simulation[END_REF][START_REF] Benioug | Numerical efficiency assessment of IB-LB method for 3D pore-scale modeling of flow and transport[END_REF].

Slip boundary condition

When LBM is applied for flow simulations in porous media where Kn > 0.001, gas slippage effect can be taken into account by adopting a slip boundary condition. Various approaches permitting to achieve this goal can be employed [START_REF] Ansumali | Kinetic boundary conditions in the lattice Boltzmann method[END_REF][START_REF] Succi | Mesoscopic modeling of slop motion at fluid-solid interfaces with heterogeneous catalysis[END_REF][START_REF] Tang | Lattice Boltzmann method for gaseous microflows using kinetic theory boundary conditions[END_REF][START_REF] Chai | Lattice Boltzmann simulation of surface roughness effect on gaseous flow in a microchannel[END_REF]. In our study, the diffuse reflection combined with the bounce-back boundary condition (DBB) is adopted following [START_REF] Chai | Lattice Boltzmann simulation of surface roughness effect on gaseous flow in a microchannel[END_REF][START_REF] Chai | Gas Flow Through Square Arrays of Circular Cylinders with Klinkenberg Effect: A Lattice Boltzmann Study[END_REF]

f i (r, t + δt) = rKf eq i (ρ w , u w ) + (1 -r) fi (r, t), if (c i -u w ) • n > 0 ( 18
)
where ρ w is the gas density at the wall, u w is the wall velocity, n is the unit wall normal, 0 ≤ r ≤ 1 is the combination coefficient, and K is the coefficient defined as

K = ξ i •n<0 |ξ i • n| fi (r, t) ξ i •n>0 |ξ i • n|f eq i (ρ w , u w ) (19)
where ξ i = c i -u w , and the incident discrete velocity c i can be expressed as

c i = c i -2(c i • n)n (20)
According to [START_REF] Chai | Gas Flow Through Square Arrays of Circular Cylinders with Klinkenberg Effect: A Lattice Boltzmann Study[END_REF], this boundary condition is more adapted to be used in porous media since it requires only local information.

For simulation of the gas slip flow, the relaxation time τ ν must be related to the Knudsen number. The Knudsen number can be related to the viscosity as [START_REF] Guo | Physical symmetry, spatial accuracy, and relaxation time of the lattice Boltzmann equation for microgas flows[END_REF]

Kn = π 2 ν c s H (21)
Then, using (8a), τ ν can be expressed as [START_REF] Guo | Physical symmetry, spatial accuracy, and relaxation time of the lattice Boltzmann equation for microgas flows[END_REF]

τ ν = 1 2 + 6 π H δx Kn ( 22
)
The relaxation time τ q and the combination coefficient r can be tuned to comply with a prescribed slip boundary condition. Usually [START_REF] Chai | Gas Flow Through Square Arrays of Circular Cylinders with Klinkenberg Effect: A Lattice Boltzmann Study[END_REF], this is done by considering the stationary Poiseuille flow between the flat walls y = 0 and y = H driven by the body force ρ(a x , 0, 0) with a slip boundary condition prescribed at the walls. In this case, the analytical solution for the parabolic velocity profile u y is given by

u y u c = 4 y H (1 - y H ) + u s u c ( 23
)
where

u c = axH 2 8ν
is the velocity at the center y = H 2 for the no-slip boundary, and u s is the slip velocity at the wall. The second order slip boundary condition can be written as [START_REF] Hadjiconstantinou | Comment on Cercignani's second-order slip coefficient[END_REF]

u s = C 1 λ ∂u ∂y wall -C 2 λ 2 ∂ 2 u ∂y 2 wall ( 24
)
where C 1 and C 2 are the coefficients describing the gas-wall interaction and can be supplied by the kinetic theory [START_REF] Hadjiconstantinou | Comment on Cercignani's second-order slip coefficient[END_REF][START_REF] Suga | Corrigendum: Lattice Boltzmann methods for complex micro flows: applicability and limitations for practical applications[END_REF] or experiments. The usual first order slip boundary condition corresponds to C 2 = 0. Following [START_REF] Hadjiconstantinou | Comment on Cercignani's second-order slip coefficient[END_REF], C 1 = 1.11 and C 2 = 0.61 values are adopted in this paper due to a good agreement with the solution of the linearized Boltzmann equation up to Kn ≈ 0.4. Using [START_REF] Guo | Analysis of lattice Boltzmann equation for microscale gas flows: Relaxation times, boundary conditions and the Knudsen layer[END_REF], the condition ( 24) can be simplified as

u s u c = 4C 1 Kn + 8C 2 Kn 2 (25)
On the other hand, by using the relation between the velocities at two fluid LBM nodes next to the wall and ( 23), u s can be expressed as [START_REF] Guo | Analysis of lattice Boltzmann equation for microscale gas flows: Relaxation times, boundary conditions and the Knudsen layer[END_REF][START_REF] Chai | Gas Flow Through Square Arrays of Circular Cylinders with Klinkenberg Effect: A Lattice Boltzmann Study[END_REF]

u s u c = 4r 2 -r 6 π Kn + 32(τ ν -1 2 )(τ q -1 2 ) -6 π(τ ν -1 2 ) 2 Kn 2 (26)
Comparing ( 25) and( 26), the combination coefficient r and the relaxation time τ q can be related to the C 1 and C 2 parameters and relaxation time τ ν as [START_REF] Chai | Gas Flow Through Square Arrays of Circular Cylinders with Klinkenberg Effect: A Lattice Boltzmann Study[END_REF]

r = 2C 1 6 π + C 1 (27a) τ q = 1 2 + 3 + π(2τ ν -1) 2 C 2 8(2τ ν -1) (27b)

Validation tests

Before tackling the complex issue of comparing with real experimental data, preliminary validations are to be conducted to assess the solution accuracy. The effect of the mesh refinement and the slip boundary condition are investigated for different geometries. In the following, all the values are given in lattice units unless the units are stated explicitly.

Poiseuille flow in a horizontal or inclined channel

The simplest way to validate the code is to consider a slip flow between two horizontal parallel walls, y = 0 and y = H, driven by a body force F = (F x , 0, 0). The analytical solution in this case is given by [START_REF] Guo | Analysis of lattice Boltzmann equation for microscale gas flows: Relaxation times, boundary conditions and the Knudsen layer[END_REF]. In order to estimate the influence of discretization, simulations are performed for various values of H ranging from 5 to 41. The (DBB) boundary condition [START_REF] Dymitrowska | Two-phase flow pore-size simulations in Opalinus clay by the lattice Boltzmann method[END_REF] is applied. The obtained dimensionless velocity profiles uy uc for Kn = 0.1 are displayed in Fig. 1. If the bounce-back boundary condition [START_REF] Humières | Multiple-relaxation-time lattice Boltzmann models in three dimensions[END_REF] combined with Eq. ( 9) results in exact velocity profile for the classical Poiseuille solution (i.e., with no-slip condition) [START_REF] Chun | Interpolated boundary condition for lattice Boltzmann simulations of flows in narrow gaps[END_REF], in contrast, a small but non-zero error persists in the presence of slippage at the solid wall. The relative error for the velocity at the center of the channel decreases from 0.018% for H = 5 to 2.93 • 10 -4 % for H = 41. Clearly, the numerical solution is very close to the analytical profile even for H = 5. The residual error is probably due to the spurious higher-order terms in the non-equilibrium distribution f i . Note that a regularization procedure was proposed in [START_REF] Suga | Corrigendum: Lattice Boltzmann methods for complex micro flows: applicability and limitations for practical applications[END_REF] to improve the solution accuracy but numerical tests that we performed have shown that it delivers larger relative errors, at least for the range of Knudsen numbers of interest (Kn < 1). These observations suggest that the introduction of slip flow boundary do not dramatically impact the LBM model sensitivity to computational mesh when the exact geometry is known. In a second more complicated case, the channel is inclined now with 45 • angle to the main lattice directions, i.e., it is aligned with y = x as displayed in Fig. 2(a). In practice, when the bounce-back [START_REF] Humières | Multiple-relaxation-time lattice Boltzmann models in three dimensions[END_REF] or the DBB (18) boundary condition is applied, the actual geometry corresponds to a staircase shape with the physical solid-fluid interface lying in the middle between the fluid and solid nodes. The LBM cartesian mesh generates a pseudo-roughness similar to the one caused by the voxel-based representation in 3D digitized images. From this test case, we investigate the influence of this loss of geometric information on the solution accuracy with and without slippage. First, the numerical solution of the fluid flow simulated in such geometry applying the bounce-back boundary condition can be compared with the analytical solution for the Poiseuille flow [START_REF] Chun | Interpolated boundary condition for lattice Boltzmann simulations of flows in narrow gaps[END_REF] with H defined as in Fig. 2(a). A very good agreement is observed as shown in Fig. 2(b), however, in contrast to the horizontal channel, the numerical solution is no longer exact [START_REF] Chun | Interpolated boundary condition for lattice Boltzmann simulations of flows in narrow gaps[END_REF]. From [START_REF] Chun | Interpolated boundary condition for lattice Boltzmann simulations of flows in narrow gaps[END_REF], the relative error in K BB when compared to the exact value K 0 due to this constant shift in the velocity profile can be estimated as -1 4H 2 . The relative error gets significantly larger when the DBB boundary condition is used to simulate the slip flow in this configuration, as it can be observed from Fig. 3(a) where the numerical solution for

H = 15 √ 2 , H = 31
K theo = K 0 (1 + 6C 1 Kn + 12C 2 Kn 2 ) ( 28
)
where

K 0 = H 2 12
is the absolute permeability of the channel. The numerically calculated permeabilities do not follow the parabolic law [START_REF] He | Theory of the lattice Boltzmann method: From the Boltzmann equation to the lattice Boltzmann equation[END_REF] ).

Originally, the relation ( 26) between the normalized slip velocity us uc and LBM parameters is obtained for the solid-fluid boundary located in the middle between the solid and fluid nodes. The condition is well defined, if we consider that the staircase geometry of the inclined channel is the true geometry. However, if we consider that this geometry corresponds to the inclined channel defined by the dashed lines in Fig. 2(a), the boundary condition [START_REF] Harrington | Gas flow in Callovo-Oxfordian claystone (COx): results from laboratory and field-scale measurements[END_REF] is not well defined for some LBM discrete velocities, i.e., the coefficients for the Kn and Kn 2 should contain extra terms due to the actual position of the solid-fluid interface as can be found in [START_REF] Tao | Boundary condition for lattice Boltzmann modeling of microscale gas flows with curved walls in the slip regime[END_REF]. Thus, considering the normalized slip velocity us uc , the source and the magnitude of the error is exactly the same for all the channels inclined at the same angle regardless the distance between the walls. Therefore, the error does not effectively decrease as H increases, and, as demonstrated in Fig. 3(b), the error grows with Kn.

Obviously, the sensitivity of the solution to the wall roughness increases as the slip effects prevail. This finding suggests that a very careful attention is required for capturing the geometric features of nanoporous materials when Kn becomes large. Another important observation is that the permeability K LBM does not converge to K BB (the permeability calculated using the bounce-back boundary condition) when Kn → 0 but to a slightly smaller value K DBB , i.e., the error introduced by the DBB condition does not disappear completely when Kn → 0.

It was mentioned that this deviation from the analytical solution can be effectively corrected by using the interpolated DBB boundary condition proposed by [START_REF] Tao | Boundary condition for lattice Boltzmann modeling of microscale gas flows with curved walls in the slip regime[END_REF] which places the solid-fluid interface exactly at the positions indicated by the dashed lines in Fig. 2(a) thus eliminating the artificial walls roughness. However, this option is not available for the geometries based on the FIB-SEM images which contain some degree of uncertainty about the actual location of the pore walls and where the solid-fluid interface 10 ).

is described by the border of the cubic voxels. As a consequence, we will not consider this numerical improvement in our model hereafter but keeping in mind this source of error for further comparison with experimental measurement.

Permeability of square arrays of circular cylinders

A more sophisticated geometry which can serve as an approximation of a porous medium is a two-dimensional square array of circular cylinders. The unit cell displayed in Fig. 4 represents the approximation of the solid cylinder inside the square region filled with fluid. A spatially periodic boundary condition is applied at the unit cell sides, and the flow inside the cell is driven by the body force F = (F x , 0, 0) which acts as a macroscopic pressure gradient ∇P . When stationary flow is obtained, the permeability can be found from the Darcy law

< v >= - K µ ∇P (30)
In [START_REF] Chai | Gas Flow Through Square Arrays of Circular Cylinders with Klinkenberg Effect: A Lattice Boltzmann Study[END_REF], the authors have performed high resolution simulations to study gas flow with Klinkenberg effect through a similar square array of solid cylinders. Their numerical domain contains multiple contiguous unit cells like the one displayed in Fig. 4 with pressure difference imposed at the opposite sides of the domain while spatially periodic boundary conditions are applied at the lateral sides. Several (n = 4 [START_REF] Chai | Gas Flow Through Square Arrays of Circular Cylinders with Klinkenberg Effect: A Lattice Boltzmann Study[END_REF]) replicated unit cells inside the domain are necessary in order to obtain grid independent results. When equilibrium is reached, the gas permeability can be calculated as [START_REF] Chai | Gas Flow Through Square Arrays of Circular Cylinders with Klinkenberg Effect: A Lattice Boltzmann Study[END_REF]

K = 2Lµ < U o > p o p 2 i -p 2 o ( 31
)
where L = l • n is the length of the domain, < U o > is the averaged velocity at the outlet, p i is the gas pressure at the inlet, and p o is the pressure at the outlet. If the fluid density variation inside the domain is small, the permeability values calculated from Eqs. [START_REF] Ismail | Gas Separation Membranes. Polymeric and Inorganic[END_REF] and [START_REF] Jones | A rapid accurate unsteady-state Klinkenberg permeameter[END_REF] are very close to each other. Thus, the use of the small unit cell (Fig. 4) with spatially periodic boundary conditions and body force as a macroscopic pressure gradient allows us to decrease the computational load without loss of precision. From these simulations, we study the mesh resolution required for predicting with accuracy the permeability of this simplified porous medium. The simulations conducted by Chai et al. [START_REF] Chai | Gas Flow Through Square Arrays of Circular Cylinders with Klinkenberg Effect: A Lattice Boltzmann Study[END_REF] are used as the reference solution. The simulations are performed with the unit cell size l = 400, and the cylinder diameter D = 361. In order to estimate the influence of the discretization, smaller unit cells are also considered with l as small as l = 25 while preserving the D l ratio. The body force F = (10 -7 , 0, 0) is small enough, so that the fluid density variations are negligible. Following [START_REF] Chai | Gas Flow Through Square Arrays of Circular Cylinders with Klinkenberg Effect: A Lattice Boltzmann Study[END_REF], the Knudsen number is defined using

D pore as H H = D pore = ε 1 -ε D ( 32
)
where ε is the porosity of the unit cell. The obtained permeabilities for various Kn values are compared with the fitting curve K = K 0 (1 + 22.52Kn + 26.87Kn 2 ) from [START_REF] Chai | Gas Flow Through Square Arrays of Circular Cylinders with Klinkenberg Effect: A Lattice Boltzmann Study[END_REF] in Fig. 5(a), and the corresponding relative deviations from this curve are displayed in Fig. 5(b).

As expected, the deviation from the [START_REF] Chai | Gas Flow Through Square Arrays of Circular Cylinders with Klinkenberg Effect: A Lattice Boltzmann Study[END_REF] fitting curve grows with degradation of discretization. It can be concluded that at least 6 lattice nodes (i.e., l = 60) spanning the pore throat are necessary in order to have the permeability value within the 10% accuracy range, thus demonstrating that the method is sufficiently robust and it can be used for simulations in coarsely discretized geometries with a reasonable precision.

Application to illite sample

The availability of extremely fine images of compacted illite sample obtained by FIB-SEM [START_REF] Gaboreau | Optimization of pore-network characterization of a compacted clay material by TEM and FIB/SEM imaging[END_REF] provides us the possibility to study its transport properties by means of numerical simulation. The obtained results can be thereafter confronted with experimental results obtained for the compacted illite sample of the same origin.

Experimental results

The images provided by [START_REF] Gaboreau | Optimization of pore-network characterization of a compacted clay material by TEM and FIB/SEM imaging[END_REF] correspond to the illite extracted from the Le Puy en Velay basin. The powdered illite was purified and compacted to reach the dry density of 1.7 g/cm 3 . The composition of the illite and the preparation processes can be found in [START_REF] Gaboreau | Optimization of pore-network characterization of a compacted clay material by TEM and FIB/SEM imaging[END_REF]. However, since the transport properties were not measured by [START_REF] Gaboreau | Optimization of pore-network characterization of a compacted clay material by TEM and FIB/SEM imaging[END_REF], a centimetre size sample was recreated from the same basic material, and the permeabilities to water and gas were experimentally measured. The measurement of permeability to gas was performed in an oedometer cell under a constant normal load of 89 MPa to ensure a dry density of 1.7 g/cm 3 . Nitrogen is injected at constant pressure and the outflow at equilibrium is measured using a soap bubble flowmeter. The permeability is then calculated by using [START_REF] Jones | A rapid accurate unsteady-state Klinkenberg permeameter[END_REF]. The permeability to water is also measured by maintaining a normal load which ensures a dry density of 1.7 g/cm 3 . First, the sample is saturated by applying water pressure on both sides with a gradient of 0.2 MPa, until the equilibrium is obtained. Then, the water is injected at constant pressure and the outflow is measured in stationary conditions to obtain the permeability by using [START_REF] Ismail | Gas Separation Membranes. Polymeric and Inorganic[END_REF].

The measured water permeability is K H2O = 2.14 • 10 -20 m 2 . The gas permeability values measured for the entry gas pressure ranging from 1.738 MPa to 12.298 MPa are displayed in Fig. 6(a) as function of 1 Pm . Using the least squares fitting, the experimental data can be approximated by ( 2 permeability starts to grow significantly which can be attributed to the modification of the initial pore space geometry. Therefore, no measurements at higher pressures where performed. The obtained K 0 value is 2.24 times higher than the K H2O value.

The differences can be attributed to the limitations of the experimental equipment and the employed measurement technique as well as uncertainties of the measured values.

Another possibility is the specific interaction of water molecules with solid walls inside the nano-pores of clayey material. Indeed, according to [START_REF] Mcphee | Core analysis: a best practice guide[END_REF], the absolute permeability for clay-rich rocks estimated from the Klinkenberg correlation ( 2) is routinely 2 -3 times higher than the permeability measured to a liquid. Another series of gas permeability measurements on a different compacted illite sample was performed for relatively low entry gas pressures (p i < 1.6 MPa). A non-linear behaviour of the measured permeability which does not follow the linear correlation (2) is observed. A similar behaviour was observed for COx argillite sample in [START_REF] Boulin | Expérimentation et Modélisation du Transfert d'hydrogène à travers des argiles de centre de stockage de déchets radioactifs[END_REF] when p i < 0.86 MPa.

Numerical results

Image sampling

Image segmentation was performed by Gaboreau et al [START_REF] Gaboreau | Optimization of pore-network characterization of a compacted clay material by TEM and FIB/SEM imaging[END_REF] using different methods: an advanced method that allows the recovery of porosity and pore distribution (see Fig. 7(a)), an erosion algorithm (see Fig. 7(b)), and the Otsu algorithm (see Fig. 7(c)). The Otsu algorithm is an often used segmentation method, in the case of the illite sample, it allows a good determination of the largest pores, essential to fluid flow [START_REF] Gaboreau | Optimization of pore-network characterization of a compacted clay material by TEM and FIB/SEM imaging[END_REF]. The erosion method favors pore connectivity and lessens isolated pores. The method favored by Gaboreau et al [START_REF] Gaboreau | Optimization of pore-network characterization of a compacted clay material by TEM and FIB/SEM imaging[END_REF], based on watersheding, shows the best fit in regards to total porosity and pore size distribution. It should be noted, however, that the fitting between image data and macroscopic measurements of porosity and pore size distribution is only fully valid if the image is representative of the sample. Gaboreau et al [START_REF] Gaboreau | Optimization of pore-network characterization of a compacted clay material by TEM and FIB/SEM imaging[END_REF] calculate the REV (Representative Elementary Volume) based on the analysis of 2D grayscale images and obtains a REV of 1000 -1200 nm. The already segmented images are 1096 × 1095 pixels large and constitute three stacks of 180 images each, meaning that the REV, as calculated by Gaboreau et al [START_REF] Gaboreau | Optimization of pore-network characterization of a compacted clay material by TEM and FIB/SEM imaging[END_REF], is not reached in the third direction. The pixel size is 5 × 5 nm. The distance between the consecutive images is also 5 nm [START_REF] Gaboreau | Optimization of pore-network characterization of a compacted clay material by TEM and FIB/SEM imaging[END_REF], therefore, the voxel size is 5 × 5 × 5 nm. The porosity ε = 0.2928 of the Fig. 7(a) sample is larger than the porosity ε = 0.1312 of the Fig. 7(b)) sample which is very close to the porosity ε = 0.1313 of the Fig. 7(c)) sample. The availability of three image stacks allows us to study the influence of the segmentation techniques on the numerically calculated transport properties.

In order to assess the role of segmentation in regards to permeability, for simulation tests, three sub-samples, hereinafter referred to as S0A, S0B, and S0C, of the same size are extracted from the three image stacks at the same location indicated by the red rectangular in Figs. 7(a), (b), and(c), respectively. The sub-samples size is 101×101×102 voxels in order to avoid high computational costs. Several solid voxels suspended inside the pore space as well as disconnected pores were removed from the sub-samples making the pore space fully connected. The total ε t and connected ε c porosities are summarized in Tab.1. All three sub-samples are percolating and can be used for fluid flow simulation. Only S0C sub-sample demonstrates significant difference between the connected and total porosities due to lower connectivity. Before the simulation, all the sub-samples are mirrored along the chosen axis in order to assure connections between the pores at the opposite sides of the sub-sample when the spatially periodic boundary condition is applied.

SOA

Considering the simulation size restriction, it can be observed that the simulated volumes are below the REV, calculated from grayscale images. On top of that, REV actually depends on the process studied, with REV for permeability being most likely larger than porosity REV [START_REF] Zhang | Pore scale study of flow in porous media: Scale dependency, REV, and statistical REV[END_REF] or grayscale images REV. In order to compensate for these limitations, a statistical study on different sub-samples is necessary. Spatial variations of geometrical and transport properties are estimated by considering 9 additional subsamples of the same size (101×101×102 voxels) extracted from each of the three original image stacks at locations displayed in Fig. 8(a). The sub-samples in these three sets are referred to as S(1-9)A, S(1-9)B, S(1-9)C following the same logic as for S0x sub-samples. The ranges of the total and connected porosities of S(1-9)x sub-samples are summarized in Tab.2. The ε t range of S(1-9)C is about two times larger than that of S(1-9)A and S(1-9)B sets which indicates a higher degree of spatial variability of geometrical properties for the images treated by the Otsu method. The sub-samples from S(1-9)x confirm a very good connectivity of three image stacks since only S6C and S8C originating from the image stack treated by the Otsu method are not percolating.

Study in continuous flow regime

Effect of sample location.

First, the Stokes flow with the no-slip bounce-back boundary condition ( 17) is simulated, and the K BB permeability of all the 30 sub-samples along the three coordinate axes was calculated. The obtained values are displayed as function of porosity in Figs. 8(b)-(d), and their variation ranges and mean values are summarized in Tab.2. The permeability values of S0x sub-samples can be found in Tab.1. It can be observed that, generally, permeability grows with porosity which is a general trend for clays [START_REF] Neuzil | How permeable are clays and shales?[END_REF], however, significant variations related to the individual pore space geometry cannot be excluded. The trend is clearly more pronounced for the sub-samples treated by the Otsu method. For S(1-9)A and S(1-9)B sub-samples, the K BB (z) values are generally higher than K BB (y) and K BB (x). This is more obvious when considering the mean values < K BB (x) >, < K BB (y) >, < K BB (z) > which for S(1-9)A sub-samples follow in the same order as the mean local pore sizes < H x l >, < H y l >, < H z l >, thus confirming the existence of a certain spatial anisotropy which is expected for a material of this kind. S(1-9)C sub-samples do not demonstrate the same tendency with < K BB (x) > being higher then < K BB (y) > and < K BB (z) >. The permeability variations are limited to one order of magnitude for S(1-9)A sub-samples, two orders of magnitude for S(1-9)B sub-samples, and three orders of magnitude for S(1-9)C sub-samples, which confirms larger spatial heterogeneity of the samples treated by the erosion and the Otsu method. Generally, the mean permeability of the samples treated by the erosion method is smaller than that of the samples treated by the Otsu method, with < K BB > of the samples treated by the original Gaboreau et al method being the largest, with the exception of < K BB (x) >. The mean permeability of the S(1-9)B samples treated by the erosion method is closer to the experimentally calculated value, with S6B sub-sample being the closest.

Effect of image resolution.

In order to estimate the influence of discretization, a refined sample is produced by replacing each voxel of S0A sub-sample with 8 smaller voxels. The permeability of the refined sample K BB (z) = 2.0352 • 10 -18 m 2 is 20.4% smaller, which indicates that the pore space discretization effects are present but they are reasonably limited.

Study in slip flow and transitional regime

Effect of the pore morphology.

To study the Klinkenberg effect, the Kn must be related to physical conditions. For an ideal Boltzmann gas, this can be done by using the following expression of Kn

Kn = k B T √ 2πd 2 P g H (33)
where k B is the Boltzmann constant, T is the temperature, and d is the gas particle kinetic diameter. Thus for H = 53.4750 nm and the nitrogen gas (d N2 = 0.364 • 10 -9 nm [START_REF] Ismail | Gas Separation Membranes. Polymeric and Inorganic[END_REF]) at T = 293.15 K and P g = 10 5 Pa, the Knudsen number can be estimated as Kn = 1.29, which corresponds to the transitional flow regime. At the nuclear waste storage conditions, much higher gas pressures can be expected notably due to the gas accumulation. Also, pressures exceeding P g = 10 6 Pa are necessary for gas to penetrate the water saturated clay, bringing the estimated Kn down to 0.129 which is still well beyond the continuous flow limits thus justifying the usage of slip-flow boundary condition for the considered illite sample.

For simulations inside a unit cell with spatially periodic boundary conditions and a small external body force, the fluid pressure variations are very small, and, therefore, P g can be considered as constant, and it can be associated with P m . On the other hand, H is essentially a local parameter and may vary significantly through the pore-space, whereas, in the presented slip-flow tests, an unique value of H was used to calculate Kn. Therefore, following [START_REF] Zhao | Study of Gas Flow Characteristics in Tight Porous Media with a Microscale Lattice Boltzmann Model[END_REF][START_REF] Lange | Pore Scale Simulation of Transport and Electrochemical Reactions in Reconstructed PEMFC Catalyst Layers[END_REF], a local pore size H l associated with each fluid lattice node was introduced for Kn calculation, in order to estimate the resulting variation of K compared to the globally defined Kn. In all the conducted slip flow simulations, the difference between the K values obtained with the globally defined H l =< H l > and with the locally variable H l is very small (less than 1%).

The H l is calculated by measuring the distances from the fluid lattice node to the solid surface following the directions of D3Q27 discrete velocities. Then, the mean value of these 13 distances is taken as the local pore size H l . The H l distribution inside the S0x and S6B sub-samples is displayed in Fig. 9(a), while the mean values are summarized in Tab.1. The peak value of S0A corresponds to the fluid lattice nodes with 50 < H l ≤ 55 nm which means that significant part of the pores are quite narrow and poorly discretized, and the mean value is < H l >= 10.6950 lattice units (or = 53.4750 nm). S0B and S6B contain a large amount of even smaller pores which results in their reduced permeability. The H l distribution of S0C sub-sample somehow differs from the others, by having a larger proportion of pores larger than 100 nm, which can be related with different porespace structure. Comparing the mean local pore sizes < H l > x , < H l > y , and < H l > z (see Tab.1) measured only along the principal axes directions, it is clear that the pores are preferentially aligned along the z-axis for all the the sub-samples followed by < H l > y , with < H l > x being the smallest.

A significant difference in pore-space structure between the sub-samples treated by the Otsu method and those treated by the Gaboreau et al [START_REF] Gaboreau | Optimization of pore-network characterization of a compacted clay material by TEM and FIB/SEM imaging[END_REF] and the erosion method can be confirmed by considering a characteristic pore size l z , which is a 2D parameter calculated for each image as

l z =
S pore l pore [START_REF] Ladd | Lattice-Boltzmann simulations of particle-fluid suspensions[END_REF] where S pore is the area of the pore, and l pore is the length of the pore border. The l z of the S0x and S6B sub-samples is displayed in Fig. 9(b). The S0A and S0B subsamples are quite similar with difference mainly due to different porosities. The S6B sub-sample extracted from a different location presents similar characteristics with S0B.

On the other hand, S0C sub-sample is quite different with l z presenting significantly larger variations and absolute values, while the S0C porosity is smaller than that of other presented sub-samples. This indicates that S0C pores have more bulbous shape compared to the narrow slit-like pores of the sub-samples treated by the Gaboreau et al [START_REF] Gaboreau | Optimization of pore-network characterization of a compacted clay material by TEM and FIB/SEM imaging[END_REF] and the erosion method thus influencing slip flow local velocities.

The Klinkenberg effect.

To study the Klinkenberg effect, the fluid flow with the slip boundary condition [START_REF] Dymitrowska | Two-phase flow pore-size simulations in Opalinus clay by the lattice Boltzmann method[END_REF] was simulated inside the S0x and S0B sub-samples along the z-axis. Kn is ranging from 10 -4 to 0.4. The Kn values can be related to the mean gas pressure inside the pores P m by means of (33) supposing the nitrogen gas at T = 293.15 K. As it was mentioned above, the simulations conducted with H =< H l > everywhere and simulations performed with site by site calculated H l provide the permeability results within 1%, therefore, H in (33) can be replaced with < H l >.

The influence of the Knudsen number on fluid velocity can be observed by considering the distribution of the scaled velocity norm u u max inside the pore space of S0A sub-sample as displayed in Fig. 10. For relatively small Kn = 0.01, Fig. 10(a) exhibits preferential pathways where the most of the flow is concentrated. As expected, the velocity is higher in the largest pores. Conversely, for large Kn = 0.4, gas velocity increases in the narrow passages due to slippage effect and fluid flow is more uniformly distributed as can be observed in Fig. 10(b). It follows, that this mechanism contributes to regulate the overall flow inside the pore space as Kn grows.

The calculated permeability for various Kn values is displayed in Fig. 6(b) where it is plotted versus the inverse of P m . The obtained K values are very well aligned and can be linearly fitted by [START_REF] Andra Dossier | Evaluation of the feasibility of a geological repository in an argillaceous formation[END_REF] with K 0 and b summarized in Tab.1. As noted in Sec.3.1, the K 0 obtained from the linear fit is smaller than K BB (z) obtained with the bounce-back boundary condition.

The Klinkenberg factor b varies significantly from sample to sample depending on the pore-space geometry. The data are highly scattered. Therefore, multiple studies were conducted in effort to correlate the b factor with other transport properties, mainly K 0 . Usually, the correlation takes the form of the power law. Thus the b factor for oil-field cores with K 0 ranging from 10 -17 m 2 to 10 -12 m 2 was correlated by [START_REF] Heid | Study of the permeability of rocks to homogeneous fluids, American Petroleum Institute, Drilling and Productions practices[END_REF] with b = 0.11K -0.39 0 [START_REF] Lange | Pore Scale Simulation of Transport and Electrochemical Reactions in Reconstructed PEMFC Catalyst Layers[END_REF] where K 0 is expressed in m 2 unit, and b is given in Pa unit. Later, the b factor for 100 core plugs with K 0 lying in the same range (from 10 -17 m 2 to 10 -12 m 2 ) was fitted by [START_REF] Jones | A rapid accurate unsteady-state Klinkenberg permeameter[END_REF] with b = 6.9K -0.36 0 [START_REF] Luo | Theory of the lattice Boltzmann method: Lattice Boltzmann models for nonideal gases[END_REF] where K 0 is expressed in mD unit, and b is given in psi unit. In [START_REF] Jones | A laboratory study of low permeability gas sands[END_REF], the b factor for low-permeable sands with K 0 varying between 10 -19 m 2 and 10 -14 m 2 was estimated as b = 0.86K -0.33 0 [START_REF] Marschall | Characterisation of Gas Transport Properties of the Opalinus Clay, a Potential Host Rock Formation for Radioactive Waste Disposal[END_REF] with K 0 and b expressed in m 2 and Pa unit, respectively. The study of the metagraywacke plug samples with K 0 ranging from 2 × 10 -22 m 2 to 4 × 10 -20 m 2 has delivered in [START_REF] Persoff | Hydrologic characterization of reservoir metagraywacke from shallow and deep levels of The Geysers vapor-dominated geothermal system, California, USA[END_REF] the b factor values lying between 1 MPa and 20 MPa. Recently, the b factor was estimated by [START_REF] Tanikawa | Klinkenberg effect for gas permeability and its comparison to water permeability for porous sedimentary rocks[END_REF] based on the data for 30 Taiwan sedimentary rocks (330 data points in total) with K 0 ranged from 10 -20 m 2 to 10 -14 m 2 delivering the following estimation b = (0.15 ± 0.06)K -0.37±0.038 0 [START_REF] Mcnamara | Use of the Boltzmann equation to simulate lattice-gas automata[END_REF] where K 0 and b are given in m 2 and Pa unit, respectively. Finally, using the fluid flow between the parallel plates as a starting point, an analytical expression for b factor was obtained in [START_REF] Woudberg | Predicting the Permeability of Very Low Porosity Sandstones[END_REF] for granular and foam models

b = 12P m 2 -σ σ Kn ( 39
)
where σ is the tangential momentum accommodation coefficient.

For further analysis, these correlations together with the experimentally obtained b factor of the compacted illite powder sample, the numerically calculated b factor of the S0x and S6B sub-samples, and the b factor experimentally measured by Boulin [START_REF] Boulin | Expérimentation et Modélisation du Transfert d'hydrogène à travers des argiles de centre de stockage de déchets radioactifs[END_REF] for a clay sample are summarized in Fig. 11. It can be noted, that two points corresponding to the experimentally measured b factor values are located in close proximity to each other although the samples are of different origin, which is in a good agreement with general correlation established between the b factor and absolute permeability K 0 . It is remarkable that the b factors calculated for the sub-samples treated by the Gaboreau et al [START_REF] Gaboreau | Optimization of pore-network characterization of a compacted clay material by TEM and FIB/SEM imaging[END_REF] and the erosion method are in a very good agreement with other correlation trends in spite of the small number of points. The K 0 value of S6B sub-sample is the closest to the measured experimentally K 0 value of the illite sample, and its b factor is 33.45% smaller than that of the illite sample, which can be considered as a good numerical estimation. On the other hand, the b factor of S0C sample treated by the Otsu method falls out of the general trend followed by S0A, S0B, and S6B sub-samples, which can be explained by the different pore space structure (see Fig. 9(b)) with domination of bulbous shaped pores. Using [START_REF] Klinkenberg | The permeability of porous media to liquid and gases, American Petroleum Institute, Drilling and Productions practices[END_REF] along with H =< H l >, the Klinkenberg factor of the selected subsamples can be also estimated with [START_REF] Mcphee | Core analysis: a best practice guide[END_REF], however, resulting estimations will not follow the power law trend.

Conclusions

In this study, we address the question of the LBM applicability to nanoporous materials using the 3D TRT-LBM model combined with the diffuse reflection boundary condition to describe the slip flow regime. Simulations conducted on several test cases highlight the accuracy of the method provided that an accurate representation of the pore-scale geometry is available. In constrast with the continuous flow regime, as predicted by the bounce-back boundary condition, slip flow simulations are highly dependent on the morphological features such as pore surface roughness. Results indicate a decrease of the permeability for rough surfaces and this behavior is consistent with the model of apparent permeability derived by Darabi et al [START_REF] Darabi | Gas flow in ultra-tight shale strata[END_REF] using the Maxwell theory. In a second part, the LBM is applied to numerically estimate the transport properties of illite from the digital images treated by three different segmentation techniques. Also, these properties are measured experimentally for the illite sample recreated from the same basic material.

The obtained results demonstrate that the LBM is a versatile method which can be successfully applied to simulation of slip and no-slip flow inside the pore space, thus permitting the calculation of the absolute permeability as well as the estimation of the Klinkenberg factor. Reliable results can be obtained even for poorly discretized pore space with relatively small number of fluid lattice nodes spanning the pore. This analysis supports LBM as a valuable tool for understanding and predicting the slip flow effects at nanometer pore-scale. However, two significant concerns should be kept in mind when applying LBM to nano-porous rock permeability calculations:

First, choice of the segmentation method applied to the digital images appears to have a significant influence on the numerical results. The absolute permeability values calculated for the same sub-sample treated by different segmentation methods can differ by three orders of magnitude. The variability of absolute permeability due to the location of the sub-sample inside the original image is limited to one order of magnitude for the sub-samples treated by the Gaboreau et al [START_REF] Gaboreau | Optimization of pore-network characterization of a compacted clay material by TEM and FIB/SEM imaging[END_REF] method, two orders of magnitude for the sub-samples treated by the erosion method, and three orders of magnitude for the sub-samples treated by the Otsu method. Segmentation methods may also modify the shape of the pores which can influence significantly the calculated Klinkenberg factor. If the pore space structure is not strongly modified by segmentation method, the calculated b factor values follow well the existing empirical trends and a good estimation of b can be obtained from a sub-sample with K 0 close to the experimental value.

Second, direct LBM calculations at the VER size which is crucial for reliable permeability assessment can be particularly difficult to achieve for nano-porous rocks due to large pore size distribution. In the present study, we circumvent this computational size issue by performing a statistical analysis on several sub-samples so that the mean and variance of the sample permeability can be calculated.

Acknowledgments

The financial support from NEEDS MIPOR is gratefully acknowledged. This work was also partially funded by the ICEEL Carnot Institute (grant MultiEC-HYDRO). We are thankful to Stephane Gaboreau who provided the digital images of the compacted illite used in this article.

Figure 1 :

 1 Figure 1: Dimensionless velocity profiles for the Poiseuille flow between horizontal parallel walls (Kn = 0.1).

√ 2 ,

 2 and H = 63

√ 2 15 √ 2 )

 2152 is compared to the analytical solution. It also follows, that the relative error ≈ 11.45% has little variation with discretization, while it strongly depends on Kn as it can be observed in Fig.3(b)where the dimensionless permeability of the 45 • inclined channel (H = is compared with the analytical solution

Figure 2 :

 2 Figure 2: (a) Schematic view of the unit cell with a 45 • inclined channel. (b) Dimensionless velocity profile for the Poiseuille flow inside the 45 • inclined channel with the bounce-back boundary condition applied at the solid-fluid interface (H = 15√ 2

Figure 3 :

 3 Figure 3: (a) Dimensionless velocity profiles for Poiseuille flow inside the 45 • inclined channel with the DBB boundary condition applied at the solid-fluid interface (Kn = 0.1). (b) Dimensionless permeability of the 45 • inclined channel as function of Kn (H = 15 √ 2

Figure 4 :

 4 Figure 4: Schematic view of the unit cell of square array of circular cylinders.

Figure 5 :

 5 Figure 5: (a) Permeability of the cubic array of circular cylinders as function of Kn for various discretizations. (b) Permeability deviation from fitting curve provided by[START_REF] Chai | Gas Flow Through Square Arrays of Circular Cylinders with Klinkenberg Effect: A Lattice Boltzmann Study[END_REF] for various discretizations.

Figure 6 :

 6 Figure 6: (a) Experimentally measured permeabilities of the compacted illite sample versus the inverse of the mean gas pressure. (b) Dimensionless numerically calculated permeabilities of S0x and S6B sub-samples versus the inverse of the mean gas pressure. Solid lines are linear fits by (2).

Figure 7 :Figure 9 :Figure 10 :Figure 11 :

 791011 Figure 7: The first image (5480 × 5475 nm) in a stack of the illite sample images treated by (a) the Gaboreau et al advanced method [19]; (b) the erosion method; (c) the Otsu method. The pores are black and the solid matrix is white. The red rectangular (L=505 nm) indicates the image size and the location of S0x sub-samples. (d) the 3D view of S0A (505 × 505 × 510 nm) sub-sample pore space.

Table 1 :

 1 Geometrical and transport characteristics of SOx, S6B sub-samples.

	SOB	SOC	S6B

-18

1.7428 • 10 -19 1.3475 • 10 -20 5.3850 • 10 -20 K 0 , m 2 2.3556 • 10 -18 1.4875 • 10 -19 1.3030 • 10 -20 4.7242 • 10

  BB (x), m 2 8.55 • 10 -19 -6.54 • 10 -18 2.97 • 10 -20 -1.63 • 10 -181.78 • 10 -20 -1.54 • 10 -17 K BB (y), m 2 9.96 • 10 -19 -5.25 • 10 -18 5.75 • 10 -20 -7.91 • 10 -19 6.95 • 10 -20 -6.08 • 10 -18 K BB (y), m 2 1.36 • 10 -18 -8.48 • 10 -18 5.38 • 10 -20 -1.24 • 10 -18 8.91 • 10 -20 -2.90 • 10 -18 < K BB (x) >, m 2 2.2807 • 10 -18 3.6648 • 10 -19 2.7376 • 10 -18 < K BB (y) >, m 2 2.5556 • 10 -18 3.0459 • 10 -19 1.5531 • 10 -18 < K BB (z) >, m 2 3.5265 • 10 -18 4.3770 • 10 -19 8.5057 • 10 -19

		S(1-9)A	S(1-9)B	S(1-9)C
	ε t	0.2539 -0.3808	0.1099 -0.2153	0.0738 -0.2858
	ε c	0.2442 -0.3779	0.1024 -0.2142	0.0274 -0.2843
	K

Table 2 :

 2 Variation ranges and mean values of geometrical and transport characteristics for S(1-9)x sub-samples.


Keywords: Multilevel converters, IGCT, IGBT, Switching losses, Conduction losses

IGCTs and IGBTs are compared in the case of a HVDC MMC. Specific figures of merit, and a current metric providing simple means to compare them, are introduced and discussed. Simulation results of a MMC model and figures of merit are shown to provide consistent results, proving that the proposed figures of merit are a very simple and fast way to select the best semiconductor switch. Furthermore, our analysis supports the growing interest in IGCTs for MMCs, as they are found to produce the lowest level of losses.

Figures-of-Merit and current metric for the comparison of IGCTs and IGBTs in Modular Multilevel Converters

Arthur Boutry 12 , Cyril Buttay 2 , Dong Dong 3 , Rolando Burgos 3 , Bruno Lefebvre 1 , Florent Morel 1 , and Colin Davidson 4 1 Introduction

The Modular Multilevel Converter (MMC) is a Voltage Source Converter (VSC) developed and used for Medium or High Voltage Direct Current (MVDC or HVDC) applications. It is a reversible, ac/dc Converter. The MMC (fig. 1a) is based on submodules (SMs), its elementary building blocks. Mainly composed of switches and a capacitor, a submodule can be seen as a small voltage source that can be inserted or not (depending on the switching sequence and the current sign along the submodules). The main type of submodule has a half-bridge topology, consisting in two switches and their associated antiparallel diodes, one capacitor and auxiliary systems.

Most submodule designs rely on IGBTs as the semiconductor switch. This is the case, for example in HVDC-MMCs produced by General Electric, ABB, Siemens [START_REF] Sharifabadi | Design, Control, and Application of Modular Multilevel Converters for HVDC Transmission Systems[END_REF] or RXPE [START_REF] Chen | Design and testing of the world's first single-level press-pack IGBT based submodule for MMC VSC HVDC applications[END_REF]. Indeed, IGBT modules offer high voltage ratings up to 6.5 kV, high current ratings, are fully controllable with little power, and can be sourced from many suppliers. On the contrary, thyristors, which are often used in other HVDC converters would require a complex circuitry to be turned-off. In theory, high voltage SiC MOSFETs would be ideal components for such application due to their low on-state voltage and low switching losses, but their cost and current ratings are currently prohibitive. Several authors [START_REF] Guédon | IGCTs in HVDC systems: Analysis and assessment of losses[END_REF][START_REF] Zhao | Practical analytical model and comprehensive comparison of power loss performance for various MMCs based on IGCT in HVDC application[END_REF][START_REF] Modeer | Loss comparison of different sub-module implementations for modular multilevel converters in HVDC applications[END_REF] have investigated the possibility of using IGCTs (a type of gate-controlled thyristor that can be turned off without parallel-type snubber) in a MMC. They demonstrated that doing so would increase the conversion efficiency of a MMC as compared to using IGBTs. Note that IGCTs require a series-type snubber (as shown on fig. 1b) to limit the switching speed and protect the associated freewheeling diode when the IGCT turns on [START_REF] Wei | Experimental evaluation of IGCT converters with reduced di/dt limiting inductance[END_REF][START_REF] Zeng | Integrated gate commutated thyristorbased modular multilevel converters: A promising solution for high-voltage dc applications[END_REF]. Even considering the additional power dissipation caused by this snubber circuit, IGCTs were found to be more efficient than IGBTs in the case of a MMC [START_REF] Guédon | IGCTs in HVDC systems: Analysis and assessment of losses[END_REF]. In [START_REF] Ladoux | On the potential of igcts in hvdc[END_REF], the cost of the IGCT-based submodules is compared to IGBT-based submodules and is found to offer lower capital cost per megawatt and lower life-cycle cost per megawatt.

However, because of the very different principles IGBTs and IGCTs operate on, they cannot be compared directly from the figures quoted in their datasheets. Therefore, there is a concern that the advantage of IGCTs over IGBT could result from a biased comparison, using devices with very different ratings. In this paper, we introduce new Figures Of Merit (FOMs), as a means for easy and accurate selection of semiconductor switches for a MMC application. The validity of these FOMs is then assessed against a dynamic MMC model: unlike an average arm model, the current waveforms, the on-state voltages and the losses are calculated at each instant. Such model provides realistic waveforms and losses. In particular, the actual switching frequency and switching instants of the semiconductor devices are an outcome of the simulation (MMCs do not operate at a constant frequency) and not an input parameter.

Comparison between IGCTs and IGBTs

General considerations and features

Table I presents an overview of the main features and parameters of of IGCTs and IGBTs. The IGCT is the result of the evolution of the thyristor technology. The main type of thyristor is the SCR (Silicon Controlled Rectifier). It is a semiconductor device made of four layers (P + N -PN + structure) which is turned-on by applying a current pulse on its gate, and turns off naturally when current stops flowing through its anode. The GTO (Gate Turn-Off) thyristor was introduced in 1980 [START_REF] Lutz | Semiconductor Power Devices[END_REF]. Compared to the SCR, it offers a controlled turn-off capability by applying a negative voltage to the gate. The GCT is an evolution of the GTO. Two improvements are the base of the GCT: a low-inductive housing-design allowing a quick, more homogeneous turn-off without filamentation (non-uniform current distribution), and a "transparent emitter buffer layer". The buffer-layer and transparent emitter technologies consist in a weakly doped n-layer located between the n-base and the p-emitter. This results in reductions in onstate voltage and turn-off energies. The low-inductive structure of the GCT allows to rapidly redirect all the anode current from the cathode to the gate at turn-off. This prevents filamentation [START_REF] Alvarez-Hidalgo | Characterisation of 3.3kV IGCTs for Medium Power Applications[END_REF] and reduces turn-off times (typical values: 10 µs for GCTs compared to 100 µs for GTOs). Another consequence is a much better immunity to dV/dt, removing the need for a turn-off snubber. The IGCT (Integrated GCT) corresponds to a GCT switch attached to its gate drive circuit. The IGCT switching frequency is limited by its driver [START_REF] Alvarez-Hidalgo | Characterisation of 3.3kV IGCTs for Medium Power Applications[END_REF], and at constant turning-off current, the higher the switching frequency the higher the power consumption of the gate driver. Furthermore, the lifetime of the gate drive is limited by that of its electrolytic capacitors, which is affected by the operating temperature.

The IGBT is a voltage-driven semiconductor developed in the early 1980s using a combination of a MOSFET and a N-P-N transistor structure [START_REF] Baliga | Fundamentals of Power Semiconductor Devices[END_REF]. Its relatively low losses and its ability to operate at high frequencies (several kilohertz) make it widely used in many fields. It requires a simple gate unit with a low power consumption. Compared to the GCT or other thyristors, which are made on a whole semiconductor wafer, IGBT modules are formed by paralleling a number of smaller semiconductor chips, making it scalable to different power ratings.

The technologies are still being improved with new structures for the IGBT (Enhanced Trench IGBTs -TSPT+ -for ABB [START_REF] Corvasce | 3300v hipak2 modules with enhanced trench (tspt+) igbts and field charge extraction diodes rated up to 1800a[END_REF], or the Injection Enhancement Gate Transistor -IEGT -for Toshiba [START_REF] Toshiba | High-power electric solutions -system catalog[END_REF]) and the IGCT (Reverse Conducting IGCT, Reverse Blocking IGCT... [START_REF] Vemulapati | Recent advancements in IGCT technologies for high power electronics applications[END_REF]), pushing their limits with higher ratings [START_REF] Corvasce | 3300v hipak2 modules with enhanced trench (tspt+) igbts and field charge extraction diodes rated up to 1800a[END_REF][START_REF] Vemulapati | Recent advancements in IGCT technologies for high power electronics applications[END_REF][START_REF] Vemulapati | An rc-igct for application at up to 5.3kv[END_REF][START_REF] Wikstrom | The 150 mm RC-IGCT: A device for the highest power requirements[END_REF], lower on-state voltage and switching losses [START_REF] Rahimo | Future trends in high power mos controlled power semiconductors[END_REF][START_REF] Vemulapati | Recent advancements in IGCT technologies for high power electronics applications[END_REF].

Gate circuit

Large, and high power consumption [START_REF] Filsecker | Comparison of 4.5-kV press-pack IGBTs and IGCTs for medium-voltage converters[END_REF] Small, low power consumption [START_REF] Filsecker | Comparison of 4.5-kV press-pack IGBTs and IGCTs for medium-voltage converters[END_REF] Switching Frequency Lower (up to 1 kHz, typical max. value) [START_REF] Alvarez-Hidalgo | Characterisation of 3.3kV IGCTs for Medium Power Applications[END_REF] Higher (up to tens of kHz) [START_REF] Toshiba | Trends in and future outlook for semiconductor devices with enhanced energy efficiency[END_REF] 2.2 Current ratings of semiconductor switches Many values are quoted in the datasheet of a semiconductor switch regarding the on-state current it can manage. These current values and those introduced in this paper are described in table II. For the IGCT, none of them can directly and simply be linked to the MMC operation, unlike for the IGBT (for which the dc-current rating on the IGBT datasheet, I dc-igbt , corresponds to the maximum current going through the semiconductor during MMC operation, I max-semi ). As a consequence, the suitability of an IGCT for a given MMC cannot be assessed directly, and no direct comparison can be done between IGBTs and IGCTs.

Thus, a new current metric (I eq-IGCT ) has to be defined to compare efficiently IGCTs and IGBTs current capabilities. This new current metric should correspond to a generic MMC case (it is not specific to one MMC implementation in particular) and has of course to be based on data available in the IGCT datasheet. The goal is to build a current metric comparable to the dc-current quoted in IGBT datasheets I dc-igbt , considering the actual current waveforms in a MMC.

The worst case for the IGCT -i.e. when the IGCT is subject to the maximum possible current -is when the submodule is inserted for an entire ac period. Considering the asymmetric operation of the MMC (because arm currents have DC and AC components), one of the IGCTs will be subject to more current than the other. These IGCTs currents have their waveforms drawn in figure 2, assuming the current of the ac output of the MMC is perfectly sinusoidal. It can be seen that this waveform is close to a 50 Hz half-sine. The adopted approach is to calculate the average current of that waveform in the general case and link it to the maximum current value of the same waveform, i.e. finding a relation between I av-semi and I max-semi . By establishing this connection between I av-semi and I max-semi , and based on I av-IGCT we can introduce I eq-igct , which is the equivalent for the IGCT of the dc current for an IGBT (I dc-igbt , quoted in the device datasheet). The waveforms of the MMC are described in eqs. (1a) and (1b), with i a , the current on the ac-side, I arm the current through the considered arm and I d the current on the dc-side.

Here, it will be considered that the arm is the one described in equation 1b, but the same reasoning can be used with any arm and any power factor. i a =I a * sin(ωt) (1a)

I arm (t) = I d 3 + I a 2 sin ωt (1b)
A MMC can be described with the following design-related ratios (eqs. (2a) and (2b)), k being the ratio between the ac and dc currents, and m being the ratio between the ac and dc voltages. The relation between these ratios and the power factor is the eq. (2c), resulting from the hypothesis that the powers from each side (ac and dc) are equal. Typical values for m is around 0.8 ± 0.1 and for cos(φ) (power factor) is around 0.85 ± 0.1 too. So the typical value of k is around 3 ± 0.8.

k = 3I a 2I d (2a) m = 2U a U d (2b) k * m * cos(φ) =2 (2c)
The calculation of the average current through the IGCT in the worst case is described in eqs.(3a) to (3d), with the integration of the positive component of the arm current, corresponding to the current in the most loaded IGCT. The arm current is positive between t 1 and t 2 , described in (3a). The average current of the IGCT is obtained in eq. (3d).

t 1 = arcsin (-1 k ) ω ;t 2 = π -arcsin (-1 k ) ω ( 3a
)
I av-semi = 1 T t 2 t 1 I arm (t)dt = 1 T t 2 t 1 [ I d 3 + I a 2 sin ωt]dt (3b)
I a = 2k 3 I d (3c)
I av-semi = I d 6π * (π -2 arcsin (- 1 k )) + kI d 3π [cos (arcsin (- 1 k ))] (3d)
This expression of I av-semi can be simplified using a linear regression, shown in eq. (4a). This regression is done on a limited range of values of k, between 2 (minimum of k in half-bridge based MMCs according to according eq. 2c) and 10 (corresponding to normal half-bridge based MMCs, 10 being chosen high for a HVDC MMC). The regression gives a R 2 value of 0.99992.

I av-semi =I d * [α * k + β] ; α = 0.103398 ; β = 0.19019 (4a)
I max-semi is obtained from equation (1b) and equation (2a) by replacing the sinus by maximal value. The function f is then introduced in equation 5b, as the ratio between I max-semi and I av-semi .

I max-semi =I d * ( 1 + k 3 ) (5a) f (k) = I max-semi I av-semi = 1 + k 3 × (α × k + β) (5b)
An equivalent current I eq-IGCT , that can be directly compared to the I dc-igbt of an IGBT, is defined in equation 6. I eq-IGCT is equal to the minimum of I mto-IGCT (maximum current at turn-off) and I av-IGCT × f (k) to ensure it does not exceed the current that can safely be interrupted, in a normal operation. In this paper, only the normal operation is taken into account, not the behaviour in case of fault.

I eq-IGCT = min(I mto-IGCT , I av-IGCT × f (k)) (6)
f is a function that is increasing with k. Here, we consider the case in which k = 2, i.e. f (k) = 2.5116, to stay in the most general case. In this case, for many IGCTs, I av-IGCT × f (k) has a higher value than I mto-igct . Two things can be concluded from that result. First, I mto-igct can be most of the times compared to I dc-igbt (even though it needs to be checked afterwards with this method). Then, most available IGCTs have very high I mto-igct (3 to 8 kA), which means they can handle much larger currents than available IGBTs (I dc-igbt = 1 to 3 kA for >3.3 kV IGBTs) in a MMC application. This is illustrated in Fig. 3 where I dc-igbt and I eq-IGCT are plotted. Finally, I mto-igct is a limit of the IGCT related to the gate circuit (how much current the gate can divert to open the IGCT) and I av-IGCT × f (k) is a thermal related limit of the semiconductor. If I eq-IGCT is equal to I mto-igct , the limit is not related to the semiconductor itself: the semiconductor could withstand more current. 2019) and [START_REF] Filsecker | Comparison of 4.5-kV press-pack IGBTs and IGCTs for medium-voltage converters[END_REF][START_REF] Chen | Stray impedance measurement and improvement of high-power IGCT gate driver units[END_REF][START_REF] Omura | Electrical and mechanical package design for 4.5kV ultra high power IEGT with 6kA turn-off capability[END_REF].

Figures of merit for conduction and switching losses

After we established a consistent current rating for IGBT and IGCTs in the previous section, we can now use it to compare the merits of different switches. These Figures of Merit (FOMs) are related to the conduction and switching losses.

For the conduction losses, a FOM that can be used is the on-state voltage. We consider in this paper the on-state voltage model as described in eq (7a) for the following reasons: it is very accurate for IGCTs over their entire SOA and acceptable for IGBTs, even though an exponential fit would be more accurate for these devices. As this on-state voltage is dependent on the current flowing through the component, we consider the MMC waveforms to calculate an average on-state voltage value. It is equal to the on-state voltage for the average current, as it is demonstrated in eqs (7a) to (7b). Thus, this approach only uses values that are available in the datasheets, which makes it interesting for at the early stages of design, for the pre-selection of devices.

V on (I) =V 0 + R.I (7a) average(V on ) =V 0 + R.average(I) = V on (I av-semi ) (7b)
The average on-state voltage is then normalised with respect to the blocking voltage of the switch (considered with a de-rating to reach a reliability of 100 Failure-In-Time -FIT -defined in particular in [START_REF] Abb | Failure rates of igbt modules due to cosmic rays[END_REF]), to allow the comparison between semiconductors with different voltage ratings. The final figure of merit is then:

FM cond-loss (IGCT ) = V on (I av-igct ) V block-100FIT (8a) FM cond-loss (IGBT ) = V on ( I dc-igbt 2.5116 ) V block-100FIT (8b)
For the switching losses, turn-on and turn-off energies are used to compare semiconductors. In datasheets, the losses figures are always quoted at maximum rated current (I swi ), which differs from one semiconductor to another and may differ from the actual current in the application. The IGBT and IGCT switching losses can be approximated by a linear function depending on the current. For the IGBT this is only true up to the dc-current rating, but this is sufficient in the case of the MMC. To compare the switching losses, the chosen figure of merit is the following: 3 FOM validation, simulation and loss study

FM swi-loss = E datasheet V block-100FIT × I swi (9)

Model Description

As it can be seen in the literature (see Tab. III), IGCTs are always found to cause less losses than IGBTs. However, the estimated amount actually varies noticeably from one publication to another because it depends on many factors: MMC ratings, control methods, operating points, not all of which are disclosed in the corresponding publications. As a consequence, it is difficult to use published results to validate our FOMs. That is why a simple MMC model has been developed, focusing on the individual behaviour of each submodule (capacitor and semiconductors). This model is used to validate the FOM approach.

The model is written as a Matlab script. The modelling assumptions are as follows: perfect sinusoidal current in the arm and perfect sinusoidal voltage across the arm for the ac components of current and voltage; arm inductor, line resistors and energy exchange between arms not modelled; only one arm considered.

The model calculates the switching instants of the semiconductors switches with the nearest level modulation control method [START_REF] Perez | Power distribution in hybrid multi-cell converter with nearest level modulation[END_REF] and the reduced switching algorithm for capacitor voltage balancing [START_REF] Tu | Reduced switching-frequency modulation and circulating current suppression for modular multilevel converters[END_REF] (a balancing capacitor algorithm based on a voltage tolerance band). Although simple, this model produces realistic waveforms for the operating voltage and current of each submodule, which allows for a more accurate estimation of the losses than average arm models. Examples of current and voltage waveforms for the upper arm, and the capacitor voltages for some (randomly chosen) submodules are displayed in figure 4. The waveforms in figure 4b show that, as it is the case in actual MMCs, the different submodules have different switching instants, different voltages at a given time, and they respect the set voltage ripple (10% in that case) around the average value of the capacitors voltages. The waveforms in figure 4a show that the semiconductor conduction losses are calculated with realistic current waveforms. This model has been used to simulate MMC operation over different length of time, and it was found that simulating operation over 5 ac-periods was sufficient to get consistent results.

Loss Study

A simulation was realised considering the following MMC: 1 GVA, 400 submodules per arm, a mean submodule voltage of 1600 V, a dc voltage of 640 kV (pole to pole). The considered devices are: IGCTs (5SHY 35L4521) and associated diodes (5SDF 20L4520), or press-pack IGBTs (5SNA 2000K450300) which include their own diodes. Using the model, losses of 5.2 MW (0.52% of the transferred power) in the case of the IGCT and 6.72 MW (0.672%) in the case of the IGBT have been calculated. This corresponds to a loss reduction of 23.6% for the IGCT case compared to the IGBT case, which is consistent with the results presented in Tab. III. [START_REF] Zhao | Practical analytical model and comprehensive comparison of power loss performance for various MMCs based on IGCT in HVDC application[END_REF] (inverter case), with iso-loss curves drawn in orange.

Analysis of the figures of merit

The figures of merit for the ABB IGCTs and Stakpak IGBTs for different voltage ratings are displayed in Fig. 5a. They are in good agreement with the corresponding switching and conduction losses calculated using the MMC model (for 4.5 kV devices only, Fig. 5b). Both approaches confirm the superiority of IGCTs over IGBTs regarding losses in general. Fig. 5a provides a comparison which is independent from the actual characteristics of the target MMC, while Fig. 5b is specific to one particular MMC. This demonstrates the interest of these figures of merit: they provide a quick, general and efficient way to compare devices even before modelling and simulating a specific MMC.

As expected, the IGCT has lower conduction losses on average, although in some cases, with the worst IGCT and the best IGBT in terms of conduction losses, the IGBT is better. This appears both in the figure of merit and in the simulated conduction losses (see figure 5). It can be explained by the use of average on-state voltage in the figure of merit -calculated as the on-state voltage at the average current of each semiconductor -and by the use of the same MMC for the simulations -the current of this MMC is the same for all the semiconductors. This implies an under-utilisation of the IGCTs, whose on-state voltage grows slowly with the current. Indeed, as it can be seen in figure 6, the on-state voltage of an IGCT grows slower than the on-state voltage of an IGBT. In fact, the IGCT is more suitable for MMCs with higher current ratings. But it has to be kept in mind that figure 6 is a particular case of one IGCT and one IGBT.

Furthermore, the IGCT is often described as having more switching losses than the IGBT: with the figures of merit (and confirmed by the simulation), it can be seen that any considered IGCT actually has lower switching losses than any IGBT. This is due to the fact that the turn-on losses are very small for the IGCT. The biggest conduction losses difference (the best IGCT and the worst IGBT in terms of conduction losses) is 0.2% while the biggest switching losses difference is 0.055%. This is due to the low switching frequency of the devices in the case of HVDC MMCs.

Conclusion

In this paper, new metrics have been developed to compare efficiently and quickly IGCTs and IGBTs in the particular case of HVDC MMCs. Figures of merit have been proposed to allow for a fair comparison in the current ratings and losses between devices. The current rating has shown that the available IGCTs have higher usable current ratings than the available IGBTs. The proposed FOMs offer a simple, analytical way to compare given IGBTs and IGCTs for a HVDC MMC application, and to easily select the best device. The different simulations have validated the FOMs and confirmed their necessity. Indeed, not all IGCTs are superior (in terms of losses) than IGBTs, and the selection must be performed on a device-per-device basis. But the FOMs confirm an advantage -in general -of IGCTs over IGBTs regarding converter efficiency, and are consistent with more complex converter-level simulations.

Fig. 1 :

 1 Fig. 1: (a) circuit diagram of the modular multilevel converter [1] and of one of its submodules (SM); (b) Snubber circuit for an IGCT-based submodule, extracted from [2].

Fig. 2 :

 2 Fig.2: Currents in the power switches, if the submodule is kept inserted (worst case).

Fig. 3 :

 3 ABB Mitsubishi

Fig. 4 :

 4 Fig.4: Waveforms of the arm and some submodules, after simulation of the same MMC as[START_REF] Zhao | Practical analytical model and comprehensive comparison of power loss performance for various MMCs based on IGCT in HVDC application[END_REF]. (p.u. means per unit and is the instantaneous voltage of the submodule divided by the rated voltage of the submodule)

Fig. 5 :

 5 Fig. 5: (a) Figures of merit applied to ABB online catalog for IGCTs and Stakpak IGBTs (2019); (b) corresponding Switching and conduction losses for 4.5 kV ABB IGCTs and Stakpak IGBTs, simulated for the same MMC as[START_REF] Zhao | Practical analytical model and comprehensive comparison of power loss performance for various MMCs based on IGCT in HVDC application[END_REF] (inverter case), with iso-loss curves drawn in orange.

Fig. 6 :

 6 Fig.6: On-state voltage of one IGCT and one IGBT[START_REF] Filsecker | Comparison of 4.5-kV press-pack IGBTs and IGCTs for medium-voltage converters[END_REF]

Table I

 I

: Global comparison of IGCTs and IGBTs (Press-Pack) Semiconductor IGCT IGBT (Press-Pack) Snubber Needed (series, to limit turn-on dI/dt) No On-state Voltage Low (around 2V at high current) High (more than 3V at high current) Turn-on Energy loss Around 2 J (2.8 kV, 2 kA) Around 10 J (2.8 kV, 2 kA) Turn-off Energy loss Around 10 J (2.8 kV, 2 kA) Around 10 J (2.8 kV, 2 kA)

Table II :

 II Description of the current-related variables mentioned in the paper.

	Current	Switch Source	Description	Symbol
	dc-current	IGBT	Datasheet	Max. dc-current that the IGBT can conduct	I dc-igbt
	Peak Current	IGBT	Datasheet	Max. peak (1ms) current that the IGBT can switch	I pk-igbt
	Max. average on-state current	IGCT	Datasheet	Based on half-sine, no real practical meaning according to [21]	I av-igct
	Max. RMS on-state current	IGCT	Datasheet	Based on half-sine, no real practical meaning according to [21]	I rms-igct
	Max. controllable turn-off current	IGCT	Datasheet	Max. anode current that can be turned-off	I mto-igct
	Switching Current	Both	Datasheet	Condition test current for the switching energy in the datasheet	I swi
	Max. av. current MMC	Both	This paper Av. current going through the switches (worst case)	I av-semi
	Max. instantaneous
	current MMC

Both This paper Max. instantaneous current seen by the switches I max-semi Equivalent current IGCT IGCT This paper Built with the datasheet currents, equivalent of I dc-igbt for the IGCT I eq-igct

Table III :

 III Loss studies results

		Total losses Total losses	Inverter or	Loss red.
		IGBT (%)	IGCT (%)	rectifier	(%)
	[7]	0.48	0.41	Average rectifier and inverter	14.6
	[5]	0.64	0.44	Rectifier	31
	[5]	0.32	0.27	Inverter	15.6
	[6]	0.76	0.74	Rectifier	2.6
	[6]	0.83	0.64	Inverter	22.9

Acknowledgements

This work was supported by a grant overseen by the French National Research Agency (ANR) as part of the "Investissements d'Avenir" Program ANE-ITE-002-01.