Joffrey Bluthé

Benoît Bary

Eric Lemarchand

Closure of parallel cracks: micromechanical estimates versus finite element computations

Keywords: Cracked media, strain concentration, crack closure, numerical simulations, homogenization theory

3D finite element simulations have been performed in order to assess the ability of four classical homogenization schemes to model the elastic behavior of solids with parallel cracks, namely the dilute, Mori-Tanaka, self-consistent and differential schemes. The cracks have been represented by right circular cylinders with aspect ratios as low as 10 -3 in the simulations whose centroids are randomly located in the REV. Special attention has been paid to the crack aspect ratio variation predicted by the different schemes, since the goal is ultimately to propose a non-linear micromechanical model of a cracked solid, taking complete crack closure into account.

The results confirm earlier studies which showed that the differential scheme was best suited for this kind of morphology when computing elastic moduli, but additionally, we show that changes in crack aperture are also accurately predicted. It is however noted that the randomness in the positions of the cracks leads to significant scatter in the magnitude of the aperture variation inside a given simulation, which suggests that the cracks do not close simultaneously. As a consequence, nonlinear numerical simulations accounting for contact between the crack lips should be performed in order to completely validate a non-linear micromechanical model in alternate tension/compression loading cases.

Introduction

Every solid body can be seen as composite at some scale, with the exception of the single crystal solid. This trivial observation explains the strong interest in the determination of effective properties of heterogeneous systems developed in the past century. Several different strategies exist to undertake this challenging task. In the present paper, Eshelby-based micromechanics arguments and numerical tools are both derived. The first one is based on Eshelby's solution to the problem of a single ellipsoidal inclusion embedded in an infinite medium (Eshelby [11]); this method has the advantage of being often completely analytical, which is interesting for parametric analysis. Its main drawback is that a fair number of simplifying assumptions usually have to be introduced in order to conduct any calculation to the end, which may limit its range of applicability. The second method relies entirely on numerical tools and is one of brute force; it consists in generating the desired microstructures and in performing numerical experiments on said structures, using some finite element software. The main advantages of the second method are that its implementation is straightforward, that much more complex and realistic geometries can be investigated, and that it gives access to local information in the multiple phases considered instead of just average quantities. However, it may require a high computing power. Even creating the geometry can prove to be very challenging for current softwares, especially in a periodic context, and the resulting mesh is not always satisfactory, particularly when dealing with singular geometries such as crack tips. That being said, the recent advances in computer technology have made this method more and more available to researchers, and it is an interesting tool to investigate the consistency of analytical estimates.

One area of particular interest for many applications is the behavior of cracked rocks. A number of researchers have applied the tools of micromechanics to determine the behavior of cracked media in order to account for a coherent coupling between their mechanical behavior and their transport properties. This approach has been implemented in Dormieux and Kondo [START_REF] Dormieux | Approche micromécanique du couplage perméabilité-endommagement[END_REF] to estimate both the effective stiffness and hydraulic conductivity of rocks in the context of the self-consistent approximation, whereby the cracks of ellipsoidal shape are endowed with a fictitious permeability derived from a Poiseuille law in order to reduce the problem to a conductivity homogenization procedure. The resulting effective hydraulic conductivity is a function of both the aspect ratio of the cracks X and their uniform radius a, which allows the authors to take into account the impact of both crack propagation and crack opening/closure on the homogenized transport properties. These ideas were further investigated by Barthélémy [START_REF] Barthélémy | Effective permeability of media with a dense network of long and micro fractures[END_REF], who considered the presence of cracks at two different scales and studied the differences between their respective mathematical treatments, and by Lemarchand et al. [START_REF] Lemarchand | Micromechanics contribution to coupled transport and mechanical properties of fractured geomaterials[END_REF] and Lemarchand et al. [START_REF] Lemarchand | Tortuosity effects in coupled advective transport and mechanical properties of fractured geomaterials[END_REF]; in a similar approach, Vu et al. [START_REF] Vu | Theoretical predicting of permeability evolution in damaged rock under compressive stress[END_REF] proposed to model the impact of crack propagation on the effective permeability; finally, Levasseur et al. [START_REF] Levasseur | A micro-macro approach of permeability evolution in rocks excavation damaged zones[END_REF] built on Dormieux and Kondo [START_REF] Dormieux | Approche micromécanique du couplage perméabilité-endommagement[END_REF] in order to predict the permeability change related to damage caused by an excavation process. Following the latter, the same kind of methodology can be applied to the opposite process, i.e. permeability reduction caused by the mechanical closure of fractures during the service life of a drift. A first attempt was made in Bluthé et al. [START_REF] Bluthé | Micromechanical modeling of the Compression of the Damaged Zone experiment in the Callovo-Oxfordian formation[END_REF], who proposed an estimate of the crack closure due to either an applied mechanical load or a hydration-induced swelling process in the context of the Mori-Tanaka approximation. Restricted to a linear context, the latter is currently being extended to the non-linear case associated with the complete crack closure. Complete closure is indeed expected because of radial compression due to the swelling of both the clay core and that of the damaged zone (see Bluthé et al. [4]).

The extension to non-linearity is entirely possible without much difficulty, since the theoretical background was laid out by Deude et al. [START_REF] Deude | Micromechanical approach to nonlinear poroelasticity: application to cracked rocks[END_REF]. In this work, the average crack closure is estimated in the context of the Mori-Tanaka scheme by using the strain concentration rule in the cracks in the limit of infinitesimal aspect ratios X 1, and the closure criterion reduces to ∆X < -X 0 , where X 0 is the initial aspect ratio of the crack and ∆X = X -X 0 . The effective tangent stiffness of the cracked medium is then computed as a function of which crack families are still open. Since the set of closed crack families is itself a function of the current stress state, the effective mechanical behavior of the medium is indeed non-linear. However, for such an approach to be valid, use of a particular micromechanical scheme requires it to predict accurately both the effective stiffness and the crack closure. Since the homogenization schemes rely on the estimation of the average strain concentration tensors in the different phases, one could be tempted to imagine that hypotheses that lead to good stiffness predictions might also lead to good strain concentration predictions. To the author's knowledge, this has not received much attention in the past, so this paper mainly aims at contributing on this aspect. Thus, the goal of this work is to identify a relevant homogenization scheme in terms of both the effective stiffness and the crack aperture changes during loading, in order to later use the selected scheme in a non-linear mechanical model at the scale of a structure, similarly to what can be found in the linear case in Bluthé et al. [START_REF] Bluthé | Micromechanical modeling of the Compression of the Damaged Zone experiment in the Callovo-Oxfordian formation[END_REF].

Let us first summarize briefly the literature on the confrontation of the standard micromechanical schemes with numerical simulations for cracked media. In some of the cases cited below, the authors studied both the case of randomly oriented cracks and that of parallel cracks, but we will restrict our discussion to the latter case since this is the context of the present study. The first attempts can be traced back to the mid 1990s when Huang et al. [START_REF] Huang | The numerical calculation of twodimensional effective moduli for microcracked solids[END_REF] implemented three of the most used micromechanical schemes, namely the dilute scheme, the differential scheme and the self-consistent scheme, for an elastic solid containing parallel cracks, and assessed their performance in comparison with 2D numerical simulations performed by applying the boundary element method (BEM) to a unit cell containing 25 cracks. For each value of the crack density parameter , as introduced in Budiansky and O'connell [START_REF] Budiansky | Elastic moduli of a cracked solid[END_REF], the authors generated 15 different crack distributions, which allowed for an estimation of the variability in their computation. In the studied range of crack densities, namely from 0.1 to 0.6, their numerical results lie between the dilute and the differential approximations. A few years later, Zhan et al. [START_REF] Zhan | Effective elastic moduli of two-dimensional solids with multiple cracks[END_REF] performed the same kind of numerical experiments but with another method, this time based on Muskelishvili's complex potentials. Again, their results seem to lie between the differential and the dilute approximation, but one of the numerical methods lies closer to the dilute estimate. This point, which may seem quite odd at first sight since the dilute estimate neglects any interaction between the cracks, will be discussed in section 2.2. Dahm and Becker [START_REF] Dahm | On the elastic and viscous properties of media containing strongly interacting in-plane cracks[END_REF] used both the finite element method (FEM) and the BEM to numerically compute the effective moduli of different crack distributions in 2D and found the differential estimate to fit the results better than the dilute estimate. Orlowsky et al. [START_REF] Orlowsky | Effects of parallel crack distributions on effective elastic properties -a numerical study[END_REF] confirmed this result by simulating elastic wave propagation, still in 2D. Shen and Li [START_REF] Shen | A numerical simulation for effective elastic moduli of plates with various distributions and sizes of cracks[END_REF] further established the good performances of the differential scheme using Kachanov's approximate method to solve their system of integral equations (Kachanov [18]). However, they pointed out that the discrepancy becomes significant for crack densities above 0.4. Saenger et al. [START_REF] Saenger | Effective elastic properties of randomly fractured soils: 3d numerical experiments[END_REF] and Saenger et al. [START_REF] Saenger | Effective elastic properties of fractured rocks: Dynamic vs. static considerations[END_REF] used the so-called rotated staggered grid technique to evaluate the dynamic response of cracked rocks in 3D, and their static response in 2D in the context of either uniform displacement or uniform stress conditions. The authors used this technique to explain the scatter on the numerical estimates of the moduli found in the literature and concluded that the differential scheme could not be set aside by any of the analyzed previous results. However, further 3D finite element estimates given by Grechka [START_REF] Grechka | Comparison of the non-interaction and differential schemes in predicting the effective elastic properties of fractured media[END_REF] tend to favor the dilute scheme with respect to the differential scheme (see discussion below). Note that for these numerical computations, the cracks were modeled as spheroids of aspect ratio 0.08. More recently, Charpin and Ehrlacher [START_REF] Charpin | Estimating the poroelastic properties of cracked materials[END_REF] proposed to evaluate the relevance of additional schemes in 2D, namely the interaction direct derivative estimate (IDD, see Zheng and Du [START_REF] Zheng | An explicit and universally applicable estimate for the effective properties of multiphase composites which accounts for inclusion distribution[END_REF]) and its modified version named full-range IDD (FRIDD), based on a comparison with finite element computations performed on periodic cells containing elliptical holes with an aspect ratio of either 0.1 or 0.01. The authors concluded that the FRIDD method was good at estimating the poromechanical properties of cracked media, but it should be noted that the presented results did not include the differential estimate. Finally, the most recent results are due to Vasylevskyi et al. [START_REF] Vasylevskyi | On micromechanical modeling of orthotropic solids with parallel cracks[END_REF], who performed 3D finite element calculations on 3D periodic cells of ellipsoids with an aspect ratio of 0.01 embedded in an orthotropic matrix. Here, the authors identified the differential scheme as giving the best results for out-of-plane moduli, and the Mori-Tanaka scheme for the in-plane ones. However, as pointed out by the authors, in the limit of infinitely flat spheroids the presence of these cavities does not impact the in-plane moduli, so that these results do not appear particularly relevant for choosing the appropriate scheme.

In this paper, we propose to investigate further the performance of five classical homogenization schemes in comparison with finite element simulations on periodic representative elementary volumes. As mentioned above, special attention was paid to variations in crack apertures upon loading. First, the main analytical results are recalled for the chosen approximations, namely the dilute, Mori-Tanaka, self-consistent and differential schemes and the Ponte-Castaneda and Willis bound. Both the compliance tensors and crack aperture variations are derived in all cases, and a few necessary comments on the range of applicability of the dilute scheme are made. Then, the numerical simulations are described: we specify the geometries and how they were generated, we explain how the boundary conditions and mesh size were chosen, and assess the quality of our results. Finally, the numerical results are presented and discussed, and we propose perspectives for future works.

Analytical estimates

General strategy

Consider a representative elementary volume (REV) Ω of a cracked solid with aligned parallel cracks. From a micromechanical point of view, cracks are usually assumed to be oblate spheroidal cavities-meaning that they are treated as spheroidal inclusions with volume fraction f = 4π 3 N a 2 c, where N is the number of cracks per unit volume of REV and a and c are the semi-major and semi-minor axes of the oblate spheroids respectively (see figure 1)-with a vanishing stiffness tensor (C c → 0), assumed isotropic for simplicity. Here, all the cavities are additionally assumed to be self-similar. The solid matrix itself has a volume fraction 1f and the stiffness tensor C m . The general strategy for the determination of the strain-based estimates of the effective stiffness and of the crack aperture variation ∆c is briefly recalled, following Dormieux and Kondo [START_REF] Dormieux | Micromechanics of Fracture and Damage[END_REF].

Uniform strain boundary conditions are here applied to the composite, i.e. ξ(z) = E • z is applied on the boundary ∂Ω of the REV, where z denotes the position vector within the REV. As long as the geometry is not singular, the strain tensor at any point z of the composite depends linearly on the loading parameter E through the fourth-rank concentration tensor, i.e. ε(z) = A(z) : E. Here and in the following developments, the colon denotes the standard double contraction of tensors. Note that the strain average rule ε = E implies that f A c + (1f )A m = I, where the following averaging operators have been introduced:

ε = 1 |Ω| ∫ Ω ε(z) dV; A α = 1 |Ω α | ∫ Ω α A(z) dV (α = c or m) (1)
and Ω c and Ω m denote the crack domain and the solid domain respectively. It then follows from the local constitutive equation of linear elasticity σ(z) = C(z) : ε(z) that the average stress is linearly related to the average strain, so that the homogenized stiffness tensor C hom may be defined as:

Σ = σ = C hom : E; C hom = C(z) : A(z) (2)
Using this definition, we obtain the relationship between C hom and A c :

C hom = C m + f (C c -C m ) : A c (3)
Introducing the crack density parameter = N a 3 (Budiansky and O'connell [START_REF] Budiansky | Elastic moduli of a cracked solid[END_REF]) and the aspect ratio of the cracks X = c a , the crack volume fraction reads f = 4π 3 X, so that (3) yields:

C hom = C m : I - 4π 3 XA c (4)
Within the considered assumptions, ( 4) is an exact definition of the homogenized stiffness. This exact definition cannot be used as such in practice, but the tools of micromechanics yield approximations through estimates of the average concentration tensor A c over the crack domain. Note that technically speaking, the previous analysis should be performed in the rate form because cracks satisfy the condition X 1 and therefore undergo large strains. However, as explained in Dormieux and Kondo [START_REF] Dormieux | Micromechanics of Fracture and Damage[END_REF], under this very condition XA c may be replaced by its limit as X → 0, so that the non-linearity disappears and the previous developments do hold.

The average increment of aspect ratio upon loading may then be computed:

c -c 0 c = ∆c c = (n ⊗ n) : ε c ( 5
)
where n is the unit normal to the midplane of the crack and c 0 is the value of c in the reference configuration. This quantity is singular as X → 0, but for non-propagating cracks (constant a), the average increment of aspect ratio is regular and is obtained using the strain concentration rule:

∆X = ∆c a = (n ⊗ n) : XA c : S hom : Σ (6)
Taking advantage of (4), it is convenient to express (6) in terms of S hom alone:

∆X = 3 4π (n ⊗ n) : S hom -S m : Σ (7)
which holds for any scheme and is especially convenient when the compliance tensor of the effective medium is known. It is straightforward to deduce from ( 7) that ∆X will be independent of for any scheme that predicts an affine dependence of S hom on .

Estimates of the concentration tensor A c now have to be proposed in order to make quantitative predictions of the homogenized stiffness. All the analytical micromechanical estimates are based on Eshelby's solution to the so-called transformation problem of an infinite elastic medium with stiffness tensor C 0 with a region subjected to a uniform eigenstrain (Eshelby [11]). According to this result, the strain in a single spheroidal cavity embedded in an infinite matrix which is subjected to uniform strain boundary conditions ξ(z) = E 0 • z at infinity is uniform in the inclusion and given by:

ε c = Q c,0 -1 : C 0 : E 0 ; Q c,0 = C 0 -C 0 : P c,0 : C 0 (8)
where the Hill tensors of the crack P c,0 and Q c,0 have been introduced. Both of these tensors are functions of the shape of the crack (superscript c) and of the stiffness tensor C 0 of the reference medium (superscript 0). Following the classical nomenclature (see Hill [START_REF] Hill | A self-consistent mechanics of composite materials[END_REF]), we will refer to this problem as the auxiliary problem; its solution is then used to estimate the average strain in the inclusion phase for the actual problem defined above on the composite. This is done by proposing adequate values for E 0 and C 0 .

Dilute scheme

In the dilute approximation, the interaction between the cracks is neglected, which means that each inclusion only "sees" the infinite matrix. This is obtained by setting

C 0 = C m and E 0 = E in (8), so that A c = (Q c,m ) -1 : C m . Denoting U c,m = lim X→0 X (Q c,m ) -1
, the following estimate for the effective compliance tensor is obtained:

S hom ∼ S dil = S m + 4π 3 U c,m (9)
This form together with (A.2) make it clear that only the out-of-plane Young's modulus and shear modulus are affected, which will be true of all the schemes. These moduli are given by:

       E hom 3 E ∼ E dil 3 E = 1 -16 3 (1 -ν 2 ) µ hom 13 µ ∼ µ dil 13 µ = 1 -16 3 1-ν 2-ν (10)
where E, µ and ν are the Young's modulus, shear modulus and Poisson's ratio of the isotropic matrix. This result is obtained by expansion to first order in , because of the assumption of dilution. Using ( 7) and ( 9), we get the average increment of aspect ratio:

∆X dil = U c,m 3333 Σ 33 (11)

Mori-Tanaka scheme

The Mori-Tanaka scheme is a first relevant candidate to take some interaction between the cracks into account. It is based on the assumption that each inclusion is embedded in the actual matrix, so that C 0 = C m , but that the composite with a single inclusion is subjected to the average strain in the matrix so that E 0 = ε m . Using the strain average rule to determine E 0 as a function of E, the following estimate of the effective compliance tensor is obtained:

S hom ∼ S MT = S m + 4π 3 U c,m (12)
This is formally equivalent to ( 9), but there is no need for an expansion to first order in this time, so that the effective moduli read:

         E hom 3 E ∼ E MT 3 E = 1 + 16 3 (1 -ν 2 ) -1 µ hom 13 µ ∼ µ MT 13 µ = 1 + 16 3 1-ν 2-ν -1 (13)
Since the compliance estimates are identical, [START_REF] Dahm | On the elastic and viscous properties of media containing strongly interacting in-plane cracks[END_REF] shows that the Mori-Tanaka estimate of average increment of aspect ratio is given by [START_REF] Eshelby | The Determination of the Elastic Field of an Ellipsoidal Inclusion, and Related Problems[END_REF].

Self-consistent scheme

It is usually considered that the Mori-Tanaka scheme can only take into account a limited amount of interaction between the inclusions. Stronger interactions can be obtained in other ways, the first of which is the self-consistent approximation. It consists in assuming that each crack is embedded in the sought homogenized medium itself, i.e. C 0 = C hom C sc , which yields an implicit equation for C sc . Different kinds of self-consistent estimates exist in the literature, but for simplicity we chose to follow Hashin [START_REF] Hashin | The differential scheme and its application to cracked materials[END_REF] and treated the composite as a dispersion in the present study, leading to the following implicit equation for the effective compliance:

S sc = S m + 4π 3 U c,sc (14)
Note the fundamental difference between ( 12) and ( 14) which lies in the superscripts, so that [START_REF] Hill | A self-consistent mechanics of composite materials[END_REF] does not give the homogenized compliance explicitly. This implicit equation is easily solved for E sc 3 and µ sc 13 numerically using a standard software. Once the homogenized stiffness is known, the average increment of aspect ratio can be computed using ( 7) and ( 14):

∆X sc = U c,sc 3333 Σ 33 (15)

Differential scheme

As mentioned in Hashin [START_REF] Hashin | The differential scheme and its application to cracked materials[END_REF], the main criticism of the self-consistent scheme is the existence of a so-called percolation threshold for certain kinds of porous media, that is, a certain finite porosity or crack density for which C sc is no longer positive definite. It is the case for spherical voids as well as randomly oriented penny-shaped cracks, and may be inaccurate since it seems to contradict several experimental results. The differential scheme has been proposed to account for a strong interaction between the cracks without percolation effects. It consists in building the cracked material incrementally by adding an infinitesimal amount of cracks to the material, until the desired crack density is reached. This process leads to the following differential equation for the effective compliance:

d de S di f = 4π 3 U c,di f (16)
This differential equation is easily solved for E di f

3 and µ di f

13 numerically using a standard software. The average increment of aspect ratio is then obtained using [START_REF] Dahm | On the elastic and viscous properties of media containing strongly interacting in-plane cracks[END_REF]. It is interesting to note that since 1 S di f -S m is the mean value of d de S di f (e) over [0, ], which is given by ( 16), we also have the following result:

∆X di f =< U c,di f 3333 > Σ 33 (17) with < • >= 1 ∫ 0 • de.

Pont-Castaneda and Willis bound

Another class of estimates may be obtained from the Hashin-Shtrikman variational principles. Details can be found in Ponte Castaneda and Willis [START_REF] Castaneda | The effect of spatial distribution on the effective behavior of composite materials and cracked media[END_REF], but the main points are first that they are functions of the spatial distribution of the inclusion phase in the composite, and secondly that these are explicit under the assumption of ellipsoidal symmetry of the distribution. The two simplest cases are obtained when the distribution is either spherical, or defined by the shape of the inclusions themselves. As is well known, the latter case reduces to the Mori-Tanaka estimate, and will not be discussed further, while the former case leads to the following estimates:

       E hom 3 E ∼ E PCW 3 E = 1 -240 (1-ν 2 ) 45+16 (7-15ν 2 ) µ hom 13 µ ∼ µ PCW 13 µ = 1 - 240 (1-ν) 45(2-ν)+32 (4-5ν) (18)
Again, the average increment of aspect ratio is obtained using [START_REF] Dahm | On the elastic and viscous properties of media containing strongly interacting in-plane cracks[END_REF]. An interesting feature of these estimates is that these are rigorous upper bounds for composites satisfying the spatial distribution requirements in a certain range for , which is why they are often referred to as the Ponte-Castaneda and Willis (PCW) bounds. For example, [START_REF] Kachanov | Elastic solids with many cracks: A simple method of analysis[END_REF] are rigorous upper bounds for the out-of-plane moduli of a composite with parallel cracks and a spherical distribution in the range 0 ≤ ≤ 3 4π . Finally, one should note that the effective Young's modulus and shear modulus vanish for = 45 128 and = 45 (2-ν) 16 (7-5ν) , so that they do not cover the whole range of values studied here numerically (see section 3).

Numerical simulations

Let us now turn to the numerical simulations that we performed in order to assess the performance of the previously described homogenization scheme. Although the goal is eventually to propose a non-linear model accounting for complete crack closure during compression of a cracked rock, the simulations performed here remain in a completely linear context. The effective out-of-plane Young's modulus and increment of crack aspect ratio were determined by applying a uniaxial tension to the elementary volumes so that the cracks tend to open instead of close, but compression would only reverse the sign because of linearity.

Preliminary investigations

First, let us describe the procedures used to generate the volume elements and to assess the validity of the numerical approach. The sample generation procedure is based on a random distribution of inclusion centroids of prescribed shape and size in a box, and has been applied to various studies and contexts, see e.g. Bary et al. [START_REF] Bary | Numerical and analytical effective elastic properties of degraded cement pastes[END_REF], Honorio et al. [START_REF] Honorio | Multiscale estimation of ageing viscoelastic properties of cementbased materials: A combined analytical and numerical approach to estimate the behaviour at early age[END_REF] and Bary et al. [START_REF] Bary | Analytical and 3d numerical analysis of the thermoviscoelastic behavior of concrete-like materials including interfaces[END_REF]. It makes use of the python tool Combs developed in the framework of the CAD platform Salome1 . In the present paper, following Charpin and Ehrlacher [START_REF] Charpin | Estimating the poroelastic properties of cracked materials[END_REF] and Vasylevskyi et al. [START_REF] Vasylevskyi | On micromechanical modeling of orthotropic solids with parallel cracks[END_REF], it was decided to introduce cracks with a shape approaching the one considered in the analytical developments, i.e. with a non-zero volume, and to further avoid a possible crack interpenetration by imposing a minimal distance (1/100 of the box edge dimension) between cracks, although overlapping is technically possible with Salome. To limit the mesh size, the number of cracks was restrained in the range 30 to 60, corresponding to a crack density parameter varying from 0.2 to 0.8. With the available numerical tools, it is possible to generate geometries with crack aspect ratios as low as 1/1000, to be compared with 1/100 in Vasylevskyi et al. [START_REF] Vasylevskyi | On micromechanical modeling of orthotropic solids with parallel cracks[END_REF]. As a consequence, the numerical strain/stress fields within the cracks are probably questionable due to elements whose shape is not well suited for FE simulations.

The meshes were generated from the geometries with automatic softwares plugged in Salome2 , and are composed of tetrahedrons for volumes and triangles for surfaces. A mesh refinement is systematically imposed on the crack surfaces by increasing the surface triangle density, to improve the quality of the simulation results. Note that a node is enforced on these upper and lower crack surfaces at the normal of each crack center, to determine the crack opening from numerical displacements precisely at this point. The meshes are all periodic, so as to be able to work with periodic boundary conditions in the finite element software Cast3M3 . Note that under these constraints, it proved difficult for the CAD software to generate the geometries using flat spheroids, so instead flat right circular cylinders were used, which made the generation more manageable for the software. The following points thus needed assessment:

• Mesh convergence • Representativeness of the elementary volumes • Impact of using right circular cylinders instead of spheroids The last point was particularly essential because the main original contribution of the present work is the comparison of the crack closure/opening during loading obtained from numerical simulations with that calculated with the homogenization schemes, and the geometry of the inclusions was likely to have an influence. Let us address it first.

In order to determine the effect of the geometry of the inclusions, 2D axisymmetric numerical simulations were first performed, since both spheroids and right circular cylinders possess this symmetry. These simulations consisted in a simple representation of Eshelby's problem of a single inclusion embedded in an isotropic linear elastic medium and subjected to a state of uniaxial tension Σ = Σe 3 ⊗ e 3 far from the inclusion (for numerical simulations, the width and height of the domain were 15 times the radius of the crack). Note that we took advantage of the symmetry of the problem with respect to the mid-plane of the crack and only considered half of the domain, with the symmetry condition on the displacement ξ 3 = 0 at the plane of symmetry. The convergence of the resulting crack opening with respect to the aspect ratio of the crack was analyzed for both geometries and compared to Eshelby's analytical result in the limit X → 0. Two different ways of quantifying the crack opening ∆X were tested: the domain occupied by the crack was meshed and given a vanishing Young's modulus with respect to that of the matrix (contrast: 2 • 10 -11 ) so as to be able to simply average the strain component ε 33 over the crack, and the displacement ξ 3 at the middle of the crack was also determined. The crack opening ∆X was then calculated from both results. As can be seen from figure 2, averaging the strain component does not yield satisfactory results since the right circular cylinder systematically underestimates the crack opening. Furthermore, convergence of the result as X → 0 is difficult, probably due to the quality of the mesh for X approaching 10 -3 , although an effort was made to decrease the mesh size as X got smaller. On the other hand, the data from ξ 3 at the middle of the crack is robust in terms of inclusion geometry and converges nicely as X → 0. These results validate the use of right circular cylinders but suggest to analyze crack opening/closure by looking at the relative displacements at the middle of the cracks lips instead of the average strain in the cracks. This means that meshing the inclusions is of no use here and can be avoided. More generally, we could consider more simple alternative representations of the cracks, e.g. by surfaces; this point will be investigated in a future work. Additionally, 3D simulations on an actual geometry containing 40 aligned cracks with X = 1 1000 ( = 0.4) were performed so as to check that these results transpose well to the 3D case. The out-ofplane elastic moduli obtained with the two geometries (spheroidal and cylindrical cracks) were also compared. Since the 3D geometries were generated randomly, they were not perfectly transversely isotropic, however, the deviation from transverse isotropy was very small (see below). Thus, the effective out-of-plane Young's modulus and shear modulus could be directly computed from two simulations with prescribed stress and periodic boundary conditions: one with Σ = Σe 3 ⊗ e 3 and one with Σ = Σ(e 1 ⊗ e 3 + e 3 ⊗ e 1 ). We compared these two effective moduli for identical geometries, i.e. the centroids of the cracks are at the same locations, so that the only difference is the shape of the cracks. We found a relative difference of 0.6% and 0.2% on the effective Young's moduli and shear moduli respectively, and a relative difference of 0.8% on the average values of the crack closure among the 40 cracks. These results allow resorting to cylinders instead of spheroids in the following investigations, the former being easier to work with. In order to assess the quality of the mesh used, we investigated the convergence of the effective elastic moduli and crack closure with the growing number of elements in the simulation. Since it is the crack tips that require a high number of elements, we ensured convergence on the elementary volume with the greatest number of cracks, namely 60 cracks ( = 0.8). As can be seen from table 1, the third column presents relative errors with respect to the fourth one of only a few percents. However, it took about an hour to perform the computation on the elementary volume with 4.32 million elements while it took about ten hours to perform the one with 6.29 million elements. The quantitative improvement was not worth the computation time, so for convenience, we performed all the other computations with a mesh size comparable to that of the simulation presented in the third column. Finally, let us tackle the representativeness of the elementary volumes. This was done for two values of the crack density, namely = 0.4 and = 0.8 and calculations were performed for four types of boundary conditions: periodic stress (PSBC), periodic strain (PEBC), uniform stress (USBC) and uniform strain (UEBC). Technically speaking, the solver implemented in Cast3M in the case of linear elasticity is based on a variational approach, namely the principle of minimum potential energy is applied, under specific constraints related to the boundary conditions. The surface tractions are directly included in the work of external forces, while the prescribed displacements are accounted for using Lagrange multipliers. Thus, the following functionals are stationary in the case of uniform boundary conditions:

         W(ξ ) -Σ : ∫ ∂Ω ξ (z) ⊗ n(z) dS (USBC) W(ξ ) - ∫ ∂Ω λ(z) • (ξ (z) -E • z) dS (UEBC) ( 19
)
where ξ is a virtual displacement field, Ω, ∂Ω and n(z) are respectively the domain, its boundary and its outward unit normal, W is the elastic energy, Σ and E are the applied stress and strain tensors respectively, and λ(z) represents the Lagrange multipliers. Of course these functionals are discretized in Cast3M, so that a linear system is eventually solved. Note that for uniform strain boundary conditions, the functional is stationary with respect to both the displacement field and the Lagrange multipliers. In the case of periodic boundary conditions, linear relationships between the displacement vectors on the different faces are applied. If l is the length of the side of the cubic simulation box, and assuming that the coordinate axes are chosen so that three faces of the domain correspond to z i = 0 with i ∈ {1, 2, 3}, the following functional is stationary in the case of periodic stress boundary conditions:

W(ξ ) -Σ : ∫ ∂Ω ξ (z) ⊗ n(z) dS - 3 i=1 ∫ z i =0 λ(z) • [ξ (z + l e i ) -ξ (z)] -[ξ (l e i ) -ξ (0)] dS (20)
where e i is the unit vector along the Oz i axis. The last term in [START_REF] Lemarchand | Micromechanics contribution to coupled transport and mechanical properties of fractured geomaterials[END_REF] ensures that the solution is of the form ξ(z) = E • z + ξ per (z), where ξ per (z) is the periodic part of the displacement, without having to specify the value of E, since the loading is here controlled via the stress tensor. On the other hand, the following functional is stationary in the case of periodic strain boundary conditions:

W(ξ ) - 3 i=1 ∫ z i =0 λ(z) • [ξ (z + l e i ) -ξ (z)] -E • z dS (21)
As in the case of uniform strain boundary conditions, the two functionals ( 20) and ( 21) are stationary with respect to both the displacement field and the Lagrange multipliers. Note that formally speaking, ( 20) and ( 21) are stationary for the exact same displacement field and Lagrange multipliers if the prescribed average stress and strain tensors are related by Σ = C per : E, where C per is the effective stiffness tensor of the periodic composite. Indeed, Σ and E are conjugate variables, and using one or the other is simply a matter of convenience. This is why the strain-based and stress-based periodic boundary conditions should systematically yield the same stiffness tensors. Nonetheless, they were both used here to ensure consistency of the procedures implemented in Cast3M, and so that the desired loading parameter (here the stress tensor) may be applied later on. On the contrary, the two functionals given in [START_REF] Kachanov | Effective elastic properties of cracked solids: Critical review of some basic concepts[END_REF] have nothing in common in general, and they do yield two different stiffness tensors. However, the discrepancy between the stress-based and strainbased effective stiffness tensors determined with uniform boundary conditions vanishes as the size of the volume element becomes much larger than the characteristic size of the heterogeneity. They may thus be used to determine whether or not the domain is large enough. For the four different kinds of boundary conditions, we computed the full effective stiffness tensor, and we first evaluated the departure from transverse isotropy. This was quantified by calculating the distance of the full tensor to the linear subspace of fourth-ranked transversely isotropic tensors with axis of isotropy along e 3 . The distance was defined according to the dot product on 6×6 matrices using the Mandel representations of the tensors. It can be computed easily since the projection of any symmetric fourth-ranked tensor (with no particular material symmetry) on this subspace is explicitly known as a function of its components, as is shown in Appendix B. Denoting || • || the norm associated to this dot product, and Π(C) the projection of C, the relative distance to the linear subspace considered is ||C-Π(C)|| ||C|| . The relative distances obtained are given in table 2. As can be seen, the full tensors are very close to being transversely isotropic. If the effective medium were far from being transversely isotropic, then we could conclude that the elementary volumes were not representative. Note that the converse is not true, which means that although these results are positive, they do not prove the quality of our elementary volumes. However, they do justify the introduction of the five elastic moduli of a transversely isotropic medium (see equation (A.1)). In order to truly quantify the quality of the elementary volumes, the out-of-plane Young's moduli and shear moduli were compared. As expected, the stiffest effective medium is obtained for uniform strain boundary conditions, while the softest is obtained for uniform stress boundary conditions (see table 3). As the size of the elementary volume increases, these lower and upper bounds should theoretically converge towards the same value when the volume is big enough to contain many inclusions. This limit is the actual effective stiffness of the medium. The values reported in table 3 show an important discrepancy between the uniform strain and uniform stress boundary conditions for a given value of , which suggests that the domain considered is either not big enough compared to the size of the inclusions, or does not contain enough of them. Since the periodic boundary conditions yield intermediate values, these are the boundary conditions that will be considered in the rest of this paper, despite the additional computation time required.

Results and discussion

We now present the results in terms of out-of-plane Young's modulus (E 3 ) and shear modulus (µ 13 ), as well as crack opening. The two elastic moduli were systematically determined with periodic boundary conditions with prescribed average stress, since only two simulations are then necessary, as described in section 3.1. Four values of were considered, namely 0.2, 0.4, 0.6 and The evolution with respect to of the effective Young's moduli and shear moduli are shown on figure 4 (symbols) and are compared with the analytical estimates (solid lines). For the numerical simulations, the symbols represent the average values of the moduli over the three simulations for each value of . The scatter could not be shown in this figure because it is smaller than the size of the symbol, which is evidence of its small amplitude. As can be seen, the differential scheme is the one that best fits the numerical simulations.

A couple of remarks deserve to be made at this point. First, concerning the fact that a number of authors found the dilute scheme to be the best fit (Zhan et al. [START_REF] Zhan | Effective elastic moduli of two-dimensional solids with multiple cracks[END_REF], Grechka [START_REF] Grechka | Comparison of the non-interaction and differential schemes in predicting the effective elastic properties of fractured media[END_REF]), although it is supposed to neglect interactions between the inclusions, it should be noted that these authors followed Kachanov [START_REF] Kachanov | Effective elastic properties of cracked solids: Critical review of some basic concepts[END_REF] and used stress-based estimates of the elastic moduli. The stress-based dilute estimate of the compliance tensor is indeed given by [START_REF] Dormieux | Approche micromécanique du couplage perméabilité-endommagement[END_REF], which directly yields:

E E hom 3 ∼ 1 + 16 3 (1 -ν 2 ) (22)
Given this result, it is tempting to then write E hom

3 ∼ E 1 + 16 3 (1 -ν 2 ) -1
, as the previously mentioned works did, without considering low crack densities. However, since the dilute estimate is only relevant to first order in , one should take the expansion up to first order in of the previous formula, which yields [START_REF] Dormieux | Micromechanics of Fracture and Damage[END_REF]. The same logic applies to the out-of-plane shear modulus. As explained in Dormieux and Kondo [START_REF] Dormieux | Micromechanics of Fracture and Damage[END_REF], one should be careful when using the dilute approximation, because the stress-based dilute estimate of the compliance tensor is not rigorously the inverse of the strain-based dilute estimate of the stiffness tensor. In fact, one can show that S dil : C dil -I is of the order of 2 , where I is the fourth-rank identity tensor. Thus, the previously mentioned authors actually used the Mori-Tanaka estimate instead of the dilute estimate.

Secondly, it is important to consider the distribution of the cracks and how it relates to the status of the PCW bound with respect to the present numerical results. In the present study, the crack centroids were positioned randomly in the elementary volumes, but because of the flatness of the cracks the statistical distribution of the crack phase is not isotropic. On the contrary, since the width of the cracks is much larger than their thickness and the centroids are equally spaced on average in all directions, the distances between the cracks are much lower in the (e 1 , e 2 ) plane than in the e 3 direction. Moreover, the geometries and the boundary conditions applied here are periodic, which means that the associated two-point correlation functions are also periodic. As a consequence, the hypotheses of ellipsoidal symmetry and of no long-range order, which are central to the analysis proposed in Ponte Castaneda and Willis [START_REF] Castaneda | The effect of spatial distribution on the effective behavior of composite materials and cracked media[END_REF], do not hold. The PCW bounds are therefore not rigorous bounds for the effective properties determined numerically in this work, which is why they are not respected in figure 4. Of course, the generated samples have a certain spatial distribution of cracks which results from the random positions of the centroids, but it is not so easily characterized, and one should keep in mind that the differential scheme is the best fit for this particular spatial distribution. Taking another distribution would most likely lead to a different result, so that it would be necessary to investigate the kind of distribution that is actually observed on cracked samples.

Figure 5 shows the crack openings when the elementary volume is subjected to a uniaxial tension of 1MPa. Again, there is very good agreement between the numerical results and the differential scheme. However, the scatter in the average values per simulation is much greater than for the elastic moduli, especially for = 0.8, and can here be seen. The scatter is indeed represented by the maximum (top of the symbol), minimum (bottom of the symbol) and average (central dot) values of the three simulations. This could have been expected since the moduli are macroscopic averages over the whole simulation cell while the crack opening is a very local information. Indeed, recall that the latter is estimated through the displacements of two opposite points of the crack lips, which may bias the results in case of strong crack-to-crack interaction. A direct consequence is that more realizations of samples should probably be generated and with more cracks so as to have a better characterization of this spread. Still, with the data already collected, we were able to look more closely at the statistical distribution of ∆X on the different simulations. Indeed, the values shown are the average values per simulation, but there is also scatter in each one of them, as can be seen from figures 6 and 7, where a smooth histogram is shown from combining all the simulations in the case = 0.2 and = 0.8 respectively, and vertical straight lines represent the predictions of the homogenization schemes. The predictions of the self-consistent scheme are not included because they are out of the plotted ranges. The two figures are quite different, since for = 0.2 the distribution is rather symmetric about its mean, while it is right-skewed for = 0.8. As was to be expected from figure 5, the prediction of the differential scheme is slightly above the average value for = 0.2 and slightly below for = 0.8, so that this scheme neither overestimates nor underestimates crack closure consistently. The spread in the distributions shown suggests that the hypothesis made in Deude et al. [START_REF] Deude | Micromechanical approach to nonlinear poroelasticity: application to cracked rocks[END_REF] that all the cracks close simultaneously in the non-linear model during compression should be further investigated on an actual non-linear simulation.

Conclusion

Deude et al. [START_REF] Deude | Micromechanical approach to nonlinear poroelasticity: application to cracked rocks[END_REF] proposed a non-linear micromechanical model of cracked media that takes into account complete crack closure. This approach was initially restricted to the Mori-Tanaka scheme, probably for convenience since its implementation only requires computation of the Hill tensor of a crack in an isotropic matrix. The non-linearity hinges on the computation of the closure of the cracks, which thus needs to be accurate enough not to overestimate/underestimate the stiffness of the material. Until now, the confrontation between the different micromechanical schemes and numerical simulations had focused exclusively on the effective stiffness that they predict, and not on the strain concentration in the cracks, so that it was difficult to assess how realistic the non-linear model of Deude et al. [START_REF] Deude | Micromechanical approach to nonlinear poroelasticity: application to cracked rocks[END_REF] was. The aim of this paper was precisely to shed light on this.

We focused here on five classical micromechanical estimates: the dilute, Mori-Tanaka, selfconsistent and differential schemes and the PCW bounds (section 2). Direct application of these schemes to cracked media shows that they differ in terms of both stiffness and strain concentration. The Mori-Tanaka and self-consistent schemes are the stiffest and softest respectively, and the dilute and Mori-Tanaka schemes have the particular feature that the crack closure does not depend on the crack density but only on the stiffness of the solid matrix and the mechanical load. It is thus tempting to think that these schemes will underestimate the actual crack closure.

Numerical simulations were performed using the finite element software Cast3M to investigate the validity of these results (section 3). Several periodic geometries were generated at different crack densities from 0.2 to 0.8, and were subjected to periodic boundary conditions. Consistency of the geometrical representation of the cracks was ensured by comparing different shapes with small aspect ratios, and the representativeness of the simulation box was investigated. The results show that both out-of-plane effective elastic moduli are very well predicted by the differential scheme, which confirms previous results of the literature (Dahm and Becker [START_REF] Dahm | On the elastic and viscous properties of media containing strongly interacting in-plane cracks[END_REF], Orlowsky et al. [START_REF] Orlowsky | Effects of parallel crack distributions on effective elastic properties -a numerical study[END_REF], Shen and Li [START_REF] Shen | A numerical simulation for effective elastic moduli of plates with various distributions and sizes of cracks[END_REF], Vasylevskyi et al. [START_REF] Vasylevskyi | On micromechanical modeling of orthotropic solids with parallel cracks[END_REF]). A new interesting result obtained in the present study is that this scheme is also very good at predicting the average crack opening. These results suggest that the differential scheme should be preferred to estimate both the out-of-plane moduli and crack closure of a cracked medium with parallel cracks over the Mori-Tanaka scheme, which means that it is also better suited for a non-linear model. However, it should be stressed that the nature of the sample generation process determines the spatial distribution of the phases. The present result thus only holds for randomly placed crack centroids, and actual spatial distribution in cracked media should be investigated in order to confront the present results to real geometries.

As noted in section 3.1, we have shown that it is more relevant to measure the crack aperture variation from displacements at the middle of the cracks than from the average strain in the volume of the cracks. One way to improve on our simulations would thus be to dispense with meshing the interior of the cracks, and even with giving the cracks a finite volume. Indeed, it would actually be easier to simply consider the cracks as perfectly flat, that is, as two initially coinciding surfaces which are free of surface tractions. In that case, one could contemplate using XFEM to perform the computations, since this technique can deal with such geometries. These improvements might allow us to deal with higher numbers of cracks in the simulations, which in turn should lead to more accurate results, given the difference observed here between uniform stress and uniform strain boundary conditions.

Finally, we have discussed the presence of a crack opening distribution within the simulations. This suggests that working with only its average value in a non-linear model might be insufficient, since there will exist a range of loads for which some cracks are closed and some are open. It would thus be interesting to perform non-linear numerical experiments whereby contact between the crack faces is accounted for, either explicitly if the cracks are given a finite volume, or through a fictitious non-linear spring between the crack surfaces if they are perfectly flat, as suggested above. This spring should have a negligible stiffness as long as the displacement jump does not exceed a specified value, and a very large one when this value is reached. Comparison with the non-linear analytical model whereby all the cracks close simultaneously would allow to conclude on the importance of this phenomenon. The difficulty may lie in the fact that it is not guaranteed that periodic boundary conditions lead to existence and uniqueness of the solution in this non-linear context. If uniform strain/stress boundary conditions are used instead, much larger unit cells have to be used, so that the upper and lower bounds given by the uniform strain and uniform stress boundary conditions be relatively close to each other, which could make the computation cost prohibitive.

Appendix A. Concentration tensor for a flat spheroid in a transversely isotropic matrix

Let B = (e 1 , e 2 , e 3 ) be an orthonormal basis of the Euclidean space. Let C 0 be the stiffness tensor of a transversely isotropic linear elastic medium with its axis of transverse isotropy parallel to e 3 . We will use the Mandel representation of second-rank and fourth-rank symmetric tensors, that is, both stress and strain tensors are represented by (6 × 1) column vectors denoted , where () T denotes the transpose of a row vector. The representation of fourth-rank symmetric tensors follows directly from the convention for second-rank symmetric tensors, so that both stiffness and compliance tensors have the same √ 2 and 2 factors in their representations. Then, following Hoenig [START_REF] Hoenig | The behavior of a flat elliptical crack in an anisotropic elastic body[END_REF], the representation of the compliance tensor S 0 reads:

[S 0 ] B = 1 E 1 -ν 1 -ν 2 0 0 0 -ν 1 1 -ν 2 0 0 0 -ν 2 -ν 2 1 H 0 0 0 0 0 0 1+ν 1 Γ 0 0 0 0 0 0 1+ν 1 Γ 0 0 0 0 0 0 1 + ν 1 (A.1)
The representation of the tensor U c,0 for a circular crack with unit normal e 3 can then be deduced from Hoenig [START_REF] Hoenig | The behavior of a flat elliptical crack in an anisotropic elastic body[END_REF]: where S and T are given by:

           S = 1-ν 2 1 2 (1 + ν 1 ) 1 Γ -ν 2 + (1 -ν 2 1 ) 1 H -ν 2 2 T = 1+ν 1 √ Γ (A.4)

2 1Figure 1 :

 21 Figure 1: Cracked solid Ω (left) and representation of the auxiliary problem (right).

Figure 2 :

 2 Figure 2: Crack opening assessment: spheroids versus cylinders and strain versus displacement. Here, Σ = 1MPa and E = 5GPa, ν = 0.2 in the matrix.

E 3 /

 3 E (USBC) 0.056 0.021 µ 13 /µ (USBC) 0.340 0.156 E 3 /E (UEBC) 0.389 0.170 µ 13 /µ (UEBC) 0.596 0.390 0.8, with respectively 30, 40, 50 and 60 inclusions. In each case three geometries were generated to investigate the scatter in the data. One such geometry is shown on figure3for each value of .

Figure 3 :

 3 Figure 3: Four periodic geometries generated.

Figure 4 :

 4 Figure 4: Effective elastic moduli as a function of crack density parameter. Five classical approximations versus numerical simulations.

Figure 5 :

 5 Figure 5: Crack opening as a function of crack density parameter for a 1MPa uniaxial tension. Five classical approximations versus numerical simulations.

Figure 6 :

 6 Figure 6: Crack opening distribution over the three simulations with = 0.2.

Figure 7 :

 7 Figure 7: Crack opening distribution over the three simulations with = 0.8.

  [α] B using their components in B with the convention [α] B = α 11 , α 22 , α 33 ,

Table 1 :

 1 Mesh convergence study. Parameters are identical to those of figure2. The relative error with respect to the last column is given in parentheses.

	Number of elements ×10 -6	0.51	1.83	4.32	6.29
	E 3 /E	0.0718 (35%) 0.0600 (13%) 0.0550 (3.6%) 0.0531
	µ 13 /µ	0.253 (13%)	0.237 (5%)	0.229 (1.8%) 0.225
	∆X × 10 3	1.03 (26%)	1.24 (11%)	1.35 (3%)	1.39

Table 2 :

 2 Relative distances from the subspace of transversely isotropic tensors for different boundary conditions and crack densities.

	0.4	0.8
	Rel. Dist. (PSBC) 1.5% 0.9%
	Rel. Dist. (PEBC) 1.5% 0.9%
	Rel. Dist. (USBC) 0.8% 0.7%
	Rel. Dist. (UEBC) 1.3% 1%

Table 3 :

 3 Quantitative evaluation of the representativeness of the elementary volumes. E 3 , µ 13 : effective out-of-plane Young's modulus and shear modulus. Parameters are identical to those of figure 2.

	0.4	0.8	0.4	0.8
	E 3 /E (PSBC) 0.181 0.047 µ 13 /µ (PSBC) 0.462 0.228
	E 3 /E (PEBC) 0.181 0.047 µ 13 /µ (PEBC) 0.462 0.228

https://www.salome-platform.org

http://www.distene.com

http://www-cast3m.cea.fr

Acknowledgment

The authors gratefully acknowledge financial support from CEA and EDF/ERMES.

Note that we have corrected the missing parentheses in the expression for S as found in Hoenig [START_REF] Hoenig | The behavior of a flat elliptical crack in an anisotropic elastic body[END_REF]. In the isotropic case, ν 1 = ν 2 = ν and H = Γ = 1, so that S = 1 -ν 2 and T = 1 + ν and:

Appendix B. Projection on the linear subspace of transversely isotropic fourth-ranked tensors

) be the basis of the linear space of transversely isotropic tensors with n as axis of transverse isotropy defined in Walpole [START_REF] Walpole | Fourth-rank tensors of the thirty-two crystal classes: multiplication tables[END_REF]. According to their definitions, in any orthonormal basis B = (e 1 , e 2 , e 3 ) such that e 3 = n, they have the following Mandel representations:

2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 (B.1) Note that E 3 and E 4 are unsymmetric. The rules of multiplication between these tensors are easily recovered using matrix multiplications on these representations. The dot product defined on (6 × 6) matrices is Tr A T B , where A and B are two such matrices, A T is the transpose of A and Tr(•) is the trace operator. This dot product is actually simply equal to A T :: B, where A and B are the two unsymmetric fourth-ranked tensors associated to A and B, and :: is the quadruple contraction of tensors defined as:

where Einstein's summation convention on repeated indices was used. Note that all the fourthranked tensors considered here possess the minor symmetries (i.e. A i j kl = A jikl = A i jl k ), so that they do have a Mandel representation, and the words "symmetric" and "unsymmetric" only refer to the major symmetry (A i j kl = A kli j ). The transpose of A is thus to be understood as the tensor whose components are A kli j . Its Mandel representation is obviously the transpose of the Mandel representation of A.

Using these rules, it is straightforward to show that W is orthogonal, but it is not orthonormal, since F :: F = G :: G = 2. An orthonormal basis of the linear subspace of transversely isotropic symmetric tensors can easily be obtained:

Since this basis is orthonormal, the orthogonal projection of any fourth-ranked symmetric tensor on the subspace of transversely isotropic tensors is obtained by taking its dot product with the five tensors of the basis. The components of the projection Π(C) in W are:

and the more classically used components follow:


Adrian Chmielewski

Jun Meng

Beien Zhu

Yi Gao

Hazar Guesmi

Hélène Prunier

Damien Alloyeau

Guillaume Wang

Catherine Louis

Laurent Delannoy

Pavel Afanasiev

Christian Ricolleau

Jaysen Nelayah
email: jaysen.nelayah@univ-paris-diderot.fr

Reshaping Dynamics of Gold Nanoparticles under H 2 and O 2 at Atmospheric Pressure

Keywords: Gold nanoparticles, environmental transmission electron microscopy, H 2 adsorption, O 2 adsorption, DFT, Multi Scale structure Reconstruction model

Despite intensive research eorts, the nature of the active sites for O 2 and H 2 adsorption/dissociation by supported gold nanoparticles (NPs) is still an unresolved issue 1 in heterogeneous catalysis. This stems from the absence of a clear picture of the evolution of the structural properties of Au NPs in the presence of these gases at near reaction conditions, i.e. at high pressures and high temperatures. We hereby report on the rst real-space observation of the equilibrium shapes of TiO 2 -supported model Au NPs under O 2 and H 2 at atmospheric pressure using window gas cell transmission electron microscopy (GCTEM). In situ GCTEM observations show instantaneous changes in the equilibrium shape of Au NPs under O 2 during cooling from 400 °C to room temperature. In comparison, no instant change in equilibrium shape is observed under H 2 environment. To interpret these experimental observations, the equilibrium shape of Au NPs under O 2 , atomic oxygen and H 2 gas environments was predicted using a multiscale structure reconstruction model. Excellent agreement between GCTEM observations and theoretical modelling under O 2 provides strong evidence for the molecular adsorption of O 2 on the Au NPs below 120 °C. Molecular adsorption takes place on specic Au facets which are identied in this work. In the case of H 2 , theoretical modelling predicts weak interactions with gold atoms which explain their high morphological stability under this gas. This work provides atomic structural information for the fundamental understanding of the O 2 and H 2 adsorption properties of Au NPs under real working conditions and also shows a new way to identify the active sites of heterogeneous nanocatalysts under reaction conditions by monitoring the structure reconstruction.

reactive. [START_REF] Hammer | Why Gold Is the Noblest of All the Metals[END_REF] The unique chemical reactivity of Au at the nanoscale has since the 1970's sparked great research interests in the application of supported Au NPs as heterogeneous catalysts in the chemical industry. [START_REF] Corma | Chemoselective hydrogenation of nitro compounds with supported gold catalysts[END_REF] The major breakthrough in gold heterogeneous catalysis occured in the 1980's following the seminal works of Haruta and collaborators on the oxidation of carbon monoxide. [START_REF] Haruta | Novel Gold Catalysts for the Oxidation of carbon-monoxide at a Temperature Far below 0-degrees-c[END_REF][START_REF] Haruta | Gold Catalysts Prepared by Coprecipitation for Low-temperature Oxidation of Hydrogen and of Carbon-monoxide[END_REF] Since these pioneering works, there have been numerous reports on the promotion of reactions at much lower temperatures or with much higher degrees of selectivity than over other metal catalysts using Au NPs deposited on metal oxide or carbon surfaces. [START_REF] Haruta | When Gold Is Not Noble: Catalysis by Nanoparticles[END_REF] These reactions include, among others, the oxidation of CO mentioned beforehand, [START_REF] Haruta | Gold Catalysts Prepared by Coprecipitation for Low-temperature Oxidation of Hydrogen and of Carbon-monoxide[END_REF]913 the selective oxidation reactions of alcohols, 1418 the Dussan (water-gas shift) reaction, [START_REF] Fu | Active nonmetallic Au and Pt species on ceria-based water-gas shift catalysts[END_REF] the synthesis of hydrogen peroxide [START_REF] Hutchings | Nanocrystalline gold and gold palladium alloy catalysts for chemical synthesis[END_REF][START_REF] Edwards | Switching O Hydrogen Peroxide Hydrogenation in the Direct Synthesis Process[END_REF] and the selective hydrogenation reactions. 2225 In these catalytic reactions, knowing the activation and reaction processes undergone by the reacting species is important to understand the active sites of supported Au NPs.

An essential step in oxidation and reduction reactions catalyzed by supported Au NPs is the adsorption and ultimately, the dissociation of O 2 and H 2 molecules, respectively, by the metal and/or support atoms. While it has been shown that the catalytic activity of supported Au NPs depends on a number of key parameters such as temperature, [START_REF] Haruta | Novel Gold Catalysts for the Oxidation of carbon-monoxide at a Temperature Far below 0-degrees-c[END_REF][START_REF] Haruta | Gold Catalysts Prepared by Coprecipitation for Low-temperature Oxidation of Hydrogen and of Carbon-monoxide[END_REF] particle size and shape [START_REF] Mohr | Identication of active sites in goldcatalyzed hydrogenation of acrolein[END_REF]2628 and the contact structure of the nanoparticles with the support, [START_REF] Arrii | Oxidation of CO on gold supported catalysts prepared by laser vaporization: Direct evidence of support contribution[END_REF][START_REF] Zanella | Characterization and reactivity in CO oxidation of gold nanoparticles supported on TiO 2 prepared by depositionprecipitation with NaOH and urea[END_REF] fundamental questions remain on oxidation reactions with O 2 as the oxidant and hydrogenation reactions with H 2 as the reducer. In particular, there is no general consensus on the nature of the active adsorption/dissociation sites and their localization in both types of reactions as theoretical and experimental studies undertaken on these subjects have come to conicting conclusions, with supports found for dierent perspectives. For instance, in the dissociation of H 2 molecules by Au NPs supported on stoichiometric and reduced TiO 2 surfaces, density functional theory calculations predict that the active sites are the low-coordinated neutral gold atoms located at corner or edge sites with no direct bonding to the support. [START_REF] Boronat | Active Sites for H 2 Adsorption and Activation in Au/TiO 2 and the Role of the Support[END_REF][START_REF] Boronat | Unravelling the Nature of Gold Surface Sites by Combining IR Spectroscopy and DFT Calculations. Implications in Catalysis[END_REF] The dissociative adsorption of hydrogen limited to Au atoms at corner and edges positions was also suggested following experimental studies of Au/Al 2 O 3 catalysts using a combination of in situ X-ray absorption spectroscopy, chemisorption, and hydrogen-deuterium exchange experiments. [START_REF] Bus | Hydrogen chemisorption on Al 2 O 3 -supported gold catalysts[END_REF] In contrast, other experimental works report that the perimeter of the interface between gold NPs and their supports are the active sites for hydrogen dissociation. [START_REF] Fujitani | Hydrogen Dissociation by Gold Clusters[END_REF][START_REF] Luza | Revealing Hydrogenation Reaction Pathways on Naked Gold Nanoparticles[END_REF] Similarly, in the interaction of supported Au NPs with O 2 , both low-coordinated Au atoms located at corner or edge positions [START_REF] Lopez | On the origin of the catalytic activity of gold nanoparticles for low-temperature CO oxidation[END_REF][START_REF] Mills | Oxygen adsorption on Au clusters and a rough Au(111) surface: The role of surface atness, electron connement, excess electrons, and band gap[END_REF][START_REF] Roldan | Critical size for O 2 dissociation by Au nanoparticles[END_REF] and metal-support sites at the perimeter 3740 are proposed as active sites for the adsorption and dissociation of O 2 molecules. The absence of a consensus on these issues is due to the lack of a clear picture of the evolution of the structural properties of Au NPs under O 2 and H 2 at near reactions conditions, i.e. at high pressures and high temperatures.

In this work, we provide the rst real-space snapshots of the morphological dynamics of titania-supported model Au NPs under O 2 and H 2 at atmospheric pressure using window gas cell environmental transmission electron microscopy (GCTEM). It has been experimentally observed and theoretically predicted that the bonding between the NP and the reactant could induce dynamic morphology changes of the NPs. 4143 For a better understanding of the reactivity of gold and gold-based nanocatalysts, a number of environmental microscopies have been used to study in situ the structural parameters under reactive gas environments over a wide range of pressure and temperature. [START_REF] Genty | Dynamic Processes on Gold-Based Catalysts Followed by Environmental Microscopies[END_REF][START_REF] Tao | Atomic-Scale Observations of Catalyst Structures under Reaction Conditions and during Catalysis[END_REF] Among these microscopies, environmental gas transmission electron microscopy has emerged as an important tool enabling real-time atomic-scale visualization of structural evolutions of nanomaterials under gas environments and temperature. 4652 It has already provided some insights into the reactivity of Au NPs with O 2 and H 2 molecules through atomic-scale and real-time visualization of particle ES under these gases. [START_REF] Giorgio | Dynamic observations of Au catalysts by environmental electron microscopy[END_REF][START_REF] Uchiyama | Systematic Morphology Changes of Gold Nanoparticles Supported on CeO 2 during CO Oxidation[END_REF][START_REF] Yoshida | Visualizing Gas Molecules Interacting with Supported Nanoparticulate Catalysts at Reaction Conditions[END_REF] In situ monitoring of the ES of supported Au NPs under O 2 and H 2 atmospheres were carried out using chemically-sensitive high-angle annular dark eld (HAADF-STEM) imaging on a JEOL JEM-ARM200F microscope operated in scanning transmission electron microscopy (STEM) mode. The gaseous media were conned in a restricted zone around the sample in the TEM using an environmental high temperature and high pressure gas cell (HTPGC) from Protochips Inc. The gas cell is a modied version of the Atmosphere TEM environmental gas cell from the same company. It incorporates two microfabricated silicon chips (so called E-chips). It can withstand static and dynamic gas environments up to atmospheric pressure and temperatures up to 1200 °C. Figure 1a shows a low magnication scanning electron microscopy image of the E-chip which serves both as particle support and heater. It is patterned with six electron transparent SiN windows for TEM observation.

Surfactant-free Au NPs were synthesized by pulsed laser deposition (PLD) in a high vacuum chamber at 10 -5 Pa. Au NPs were deposited on rutile titania (r-TiO 2 ) nanorods synthesized by hydrothermal method and exposing (110) facets. [START_REF] Li | On the selective growth of titania polymorphs in acidic aqueous medium[END_REF] TiO 2 is a good metal oxide support for Au NPs due to its strong interaction with the latter and its chemical stability among others. [START_REF] Bagheri | Titanium dioxide as a catalyst support in heterogeneous catalysis[END_REF] The (110) facet is the most stable facet of r-TiO 2 and is widely used as substrate in heterogeneous catalysis. [START_REF] Diebold | The surface science of titanium dioxide[END_REF] The nanorods were deposited prior to metal deposition on the smaller heater E-chip of the HTPGC. The low magnication (×1M) HAADF-STEM image in Figure 1b shows r-TiO 2 nanorods randomly dispersed on one of the SiN windows. The nanorods were heated at 300 °C during metal deposition and the nominal thickness of metal deposited xed to 1 nm to obtain NPs with a mean particle in-plane size around 5 nm, a size range where TiO 2 -supported Au NPs in fact start to be active. [START_REF] Fujitani | Hydrogen Dissociation by Gold Clusters[END_REF][START_REF] Haruta | Lowtemperature Oxidation of Co Over Gold Supported on TiO 2 , α -F e 2 O 2 , and Co 3 O 4[END_REF] Further experimental details on electron microscopy, the gas cell and sample preparation are provided in the Supporting Information. Besides O 2 and H 2 , Au NPs were also studied in situ under residual argon gas in a benchmark experiment. Under all gases, the Au NPs were rst heated in situ to 400 °C and their ES were imaged at dierent temperature steps while cooling the samples to RT. Figure 1c shows a higher magnication (×5M) HAADF-STEM image of an assembly of Au NPs on the r-TiO 2 support under 40 Pa Ar after in situ heating to 400 °C. The in-plane particle size measured from the HAADF-STEM image is below 6 nm. By measuring the projection angles between the edges of the projected outlines of individual NPs and comparing them with the angles between crystal planes of the Au fcc structure with dierent Miller Indices (see Figure S1 in supplementary information), we deduce that the particle outlines are bounded mainly by {111} and {100} edge-on facets pointing to a nearly truncated octahedron ES as observed in similar Au NPs under vacuum. [START_REF] Tai Nguyen | Structural Properties of Catalytically Active Bimetallic Gold-Palladium Nanoparticles Synthesized on Rutile Titania Nanorods by Pulsed Laser Deposition[END_REF] In Figure 1c, the majority of the truncated octahedra are viewed close to the [101] zone axis (like the ones pointed by white arrows). Along this direction, their STEM projection images display a 3-fold symmetry. The anisotropy ratio γ 100 /γ 111 between the surface energies of the {100} and {111} facets, determined by applying Wul construction to the projected HAADF-STEM images of 5 NPs in [101] zone axis, is 1.10 ± 0.04. In comparison, the γ 100 /γ 111 ratio determined in vacuum for TiO 2 -supported Au truncated octahedra is 1.08 ± 0.04 (Figure S2). Thus, at 400 °C under Ar gas, the Au NPs adopt a vacuum Wul-Kaishev ES.

To establish the temperature dependence of the ES, we have studied an assembly of TiO 2supported Au NPs cooled under 40 Pa Ar in a stepwise manner from 400 °C to RT with 100 °C temperature increments at a cooling rate of 1 °C/s between two consecutive steps. This cooling rate is applied in all cooling experiments presented afterhand in this work.

Figure 1d shows the corresponding temperature series of HAADF-STEM images of the NPs.

Here, the NPs are supported on both r-TiO 2 nanorods and anatase TiO 2 NPs. The anatase nanostructures are produced as residual by-products during the synthesis of rutile nanorods.

Analysis of the projected HAADF-STEM images of individual NPs on r-TiO 2 shows NPs oriented close to the [101] zone axis as well as in other zone axes along which their projected images show no particular symmetry. During the cooling process, we rst observe that, due to the strong interactions between the TiO 2 supports and the Au NPs, the latter do not diuse on the supports leading to a stable size distribution with temperature. Secondly, the NPs keep their truncated octahedron ES down to RT. This benchmark experiment shows that in the absence of an interaction between the Au NPs and its gas environment, the ES, below 400 °C, is of Wul-Kaishev type with γ 100 /γ 111 ratio equal to about 1.10. ES is driven solely by the sample temperature and particle-support interaction. Following this observation, we have performed in situ GCTEM studies of the evolution of ES of the Au NPs under O 2 and H 2 at atmospheric pressure and below 400 °C to unveil the interaction of the Au NPs with these molecules. In order to gain more detailed information of the observed morphological transformation in Figures 2c-d viewed close to [101] zone axis between 200 to 45 °C using HAADF-STEM imaging at ×10M magnication, as shown in Figures 2f-l. For comparison, the HAADF-STEM image of the NP at 300 °C where it displays a perfect truncated octahedron ES is shown in Figure 2e.

The particle outline at this temperature (yellow dotted lines on Figure 2e) and the Miller indices of the six edge-on facets along its perimeter are reproduced on all subsequent images in the temperature series to highlight deviations from the truncated octahedron ES during cooling. From the temperature series, we observe that the truncated octahedron ES is stable down to 145 °C (Figures 2e-h). At 120 °C, the NP outline starts to develop round edges (Figure 2i). The locations of these rounded smooth edges are indicated by white arrows and their extension delimited by red dotted lines on the corresponding STEM image. These rounded edges result from the extension of existing facets and/or stabilization of new ones between the adjacent {111} facets and between the {111} and {100} facets. These additional facets most probably correspond to the {110} facets and high-index facets such as {113} that connect {111} and {100} facets in a fcc structure. [START_REF] Métois | Sem Studies of Equilibrium Forms -Roughening Transition and Surface Melting of Indium and Lead Crystals[END_REF] Such facets possess a high density of low-coordinated atoms forming steps, edges and kinks. [START_REF] Quan | High-Index Faceted Noble Metal Nanocrystals[END_REF] It has been shown that these defects are catalytically active sites. [START_REF] Fujita | Atomic origins of the high catalytic activity of nanoporous gold[END_REF] In comparison, the {111} and {100} facets are ordered and at. From 120 to 45 °C, the extension of these new facets increases at the expense of the {111} and {100} ones leading to a rounded ES at 45 °C (Figures 2j-l). We also note that while particle shape changes, there is no signicant variation in the size of the particle in projection.

Figure 2m shows the variation as a function of temperature of the percentage of the particle outlines in Figures 2e-l which is delimited by either {111}-(•) or {100}-type edge-on facets (•). Between 300 and 145 °C, as the ES is stable and close to a truncated octahedron, the extension of {111}-and {100}-type edge-on facets are thus constant (with the percentage of the particle outline covered by the {111} facets being equal to 67 %). As the temperature is reduced from 145 to 45 °C and the particle ES evolves with the stabilization of new facets, we observe a decrease in the extension of both the {111}-and {100}-type facets with temperature. The rate at which the extension of the {111} facets decreases is higher than that of the {100} ones. At 45 °C, both types of edges have similar extension (equal to about 25 %). The total fraction of {111}-and {100}-type facets (•) indicates that the extension of the {111}-and {100}-type facets and the newly-developed facets are comparable at 45 °C which results in a near-perfect rounded ES in the projection at that temperature. In order to gain atomic level insights into the reshaping of Au NPs under O 2 and their high stability under H 2 , we applied the multiscale structure reconstruction (MSR) model that takes into account the direct adsorption of these gases on Au NPs facets to predict the resulting equilibrium shape. This model is based on the Wul theorem, [START_REF] Wul | On the question of speed of growth and dissolution of crystal surfaces[END_REF] the Langmuir adsorption isotherm [START_REF] Langmuir | The Adsorption of Gases on Plane Surfaces of Glass, Mica and Platinum[END_REF] and the density functional theory (DFT) calculations.

Under gas environments, and as detailed in previous works, [START_REF] Zhu | Shape Evolution of Metal Nanoparticles in Water Vapor Environment[END_REF]66 the surface tension γ hkl is corrected to be the interface tension γ int hkl by adding the contribution of adsorbents on each (hkl) surface according to:

γ int hkl = γ hkl + θ hkl E ads A at (1)
where θ hkl is the coverage of the gas on the (hkl) facet, E ads is the adsorption energy of gas molecules on this facet, and A at is the corresponding area per gold surface atom.

In this equation, the coverages θ hkl depend on the temperature T , the pressure P , and the value of E ads . We considered the Langmuir adsorption isotherm to study the O 2 @Au and the H 2 @Au systems such that:

θ hkl = P K/(1 + P K) (2)
where K is the equilibrium constant, which can be calculated as follow:

K = exp(-∆G/RT ) = exp(-E ads -(S ads -S gas )/RT ) (3)
where R is the ideal gas constant, and S gas is the gas-phase entropy of the adsorbent, obtained from the NIST-JANAF thermochemical tables. 67 The S ads represents the entropy in adsorbed state, which is assumed as zero. Based on the Wul theorem, knowing the (T , P ) dependent Based on DFT calculated interfacial energies, we modeled Au NPs with radius of 5 nm size, close to the experimental samples (around 25000 atoms) and under gas environment and temperature conditions. As shown in Figure 4, the ESs of Au NPs exposed to 1 bar of molecular oxygen are predicted to change upon cooling. For temperature higher than 17 °C, the ES of Au NP is found to be a perfect truncated octahedron as in the case of the experimental observations under vacuum and under the inert Ar gas. For temperature higher than 47 °C, the ES of Au NP is found to be a perfect truncated octahedron as in the case of the experimental observations under vacuum and under the inert Ar gas. This is because at high temperature, the coverage of O 2 on all Au facets is zero, as shown in Figure 5a.

Thus, the surface tension correction is zero and there is no reshaping. Interestingly, when the temperature decreases, the shape of Au NPs becomes more and more round. Indeed, the adsorption coverage curves shows that at low temperature, O 2 molecules adsorb on the (110) and (100) facets but not on (111) facets (see adsorption energies given in Table S2). The coverage on (110) facet is much larger than the one on (100) facet at the same temperature due to the stronger bonding of O 2 to (110) facet. The fundamental mechanism of the reshaping is that the bonding between the O 2 and the Au NPs reduces the surface tension of the latter. This reduction is stronger on (110) facets than on (100) and (111) facets due to the stronger interaction of O 2 molecule on (110) facet. Consequently, the γ int 110 becomes lower than γ int 111 and γ int 100 (Figure 5b). The resulting increased in the fraction of the (110) facet leads to the morphology of NP being more spherical (Figure 5c). Similar rounding behavior can also be found in the previous reports about the structure reconstruction of Cu and Pt NPs in water vapor environments and during CO oxidations, respectively. [START_REF] Hansen | Atom-resolved imaging of dynamic shape changes in supported copper nanocrystals[END_REF][START_REF] Vendelbo | Visualization of oscillatory behaviour of Pt nanoparticles catalysing CO oxidation[END_REF] Comparing Fig 3 and Fig 4, perfect matches between theoretical modeling and experiments are achieved. This agreement provides strong evidence that the O 2 molecules are molecularly adsorbed on the Au NPs. Note that the predicted refacetting temperature is in the range from 47 °C to 17 °C which is lower than the experimental temperatures. Although we have adopted the dispersion correction in the above calculations, the DFT-GGA method tends to underestimate energies (surface energies and adsorption energies), which is a known artifact. [START_REF] Gautier | Molecular adsorption at Pt(111). How accurate are DFT functionals[END_REF] This results in the predicted theoretical temperatures being generally lower than the experimentally measured ones. Since the dissociative adsorption of O 2 by gold cannot be disregarded in our experiment, we have also tried to interpret our experimental observations by considering the adsorption of atomic oxygen on the Au(100), Au(111) and Au(110) facets.

DFT calculations show that the adsorption of atomic oxygen is quite similar on the three low-index facets (Table S4). In the latter case, gas absorption hardly occurs on the Au NP itself. These ndings constitute key elements for the rational design of more active Au nanocatalysts for both oxidation and hydrogenation catalytic reactions. Our work provides a fundamental understanding of the structure reconstruction of TiO 2 -supported Au NPs under reactive media, which oers useful guidance for the experiments in the various related research elds.

Table S3. Adsorption energies (eV) of hydrogen molecule on Au(100), (110), and (111) surfaces.

Table S4. Adsorption energies (eV) of atomic oxygen on Au(100), (110), and (111) surfaces.

  By combining real-time in situ environmental TEM observations of particle equilibrium shape (ES) of Au NPs under gas atmospheres and multi-scale structure reconstruction (MSR) modeling, we additionally provide valuable atomic level insights on Au surface reconstruction induced by the adsorptions of O 2 and H 2 molecules. The dynamic shape changes under O 2 observed experimentally and predicted by MSR identify that the Au low-coordinated atoms can instantly and directly bind O 2 molecules at room temperatures. Under H 2 , no instant shape change is observed, indicating a weak interaction between Au NPs and H 2 molecule.

Figure 1 :

 1 Figure 1: Sample preparation for GCTEM and in situ GCTEM observation of Au NPs supported on r-TiO 2 nanorods exposing (110) surfaces under inert Ar atmosphere. (a) SEM image of the E-chip used both as sample support and heater in the Protochips gas cell. The scale bar indicates 10 µm. There are six 9 µm large windows overlaid with an electron transparent SiN lm. TEM observations are carried out at these windows. (b) HAADF-STEM image showing a random deposition of r-TiO 2 nanorods on one of the SiN windows. (c) HAADF-STEM image of Au NPs deposited by pulsed laser deposition on the r-TiO 2 nanorods and exposed to 40 Pa of Ar at 400 °C. Indexation of the projected outlines of the NPs in the HAADF-STEM image indicates that they have a truncated ES bounded by major {111} and {100} facets of the Au fcc structure. Some truncated octahedra imaged closed to the [101] zone axis are indicated by white arrows. (d) Temperature series of HAADF-STEM images of another assembly of TiO 2 -supported Au NPs cooled from 400 °C to room temperature under 40 Pa Ar. The particle size distribution and the truncated octahedron ES of the NPs are stable with temperature.

Figures

  Figures 2a-d show a temperature series of HAADF-STEM images of Au NPs on r-TiO 2 nanorods acquired just after the temperature has been set successively at 400 °C, 200 °C, 100 °C and RT under 10 5 Pa O 2 . Similar to in situ TEM observations under Ar, the NPs display truncated octahedron ES in dierent zone axes at 400 °C (Figures 2a).During sample

  During sample cooling, truncated octahedron ES is stable down to 200 °C (comparison of Figures 2a and b). The anisotropy ratio γ 100 /γ 111 , measured on ve truncated octahedra in [101] zone axis, is equal to 1.09 ± 0.04. This clearly shows that, above 200 °C, 10 5 Pa O 2 has no eect on the ES of the particles. However, when the temperature is reduced to 100 °C or RT, we observe a rounding of particle ES (Figures 2c-d) which does not occur under Ar conditions. Particle rounding is clearly visible on the projections of those NPs initially oriented closed to the [101] zone axis which lose their 3-fold symmetry in projection. The change in particle morphology with temperature is a clear indication that the surface of the Au NPs is binding O 2 molecules. It should be noted that at each temperature step, particle ES is stabilized instantaneously and it does not evolve until the temperature is reduced to the following step. Moreover, a similar rounding of particle shape is observed upon cooling under O 2 at lower pressure (133 Pa O 2 , see Figure S2 in supplementary information).

  , we have followed the change in the projected ES of individual Au NPs

Figure 2 :

 2 Figure 2: In situ HAADF-STEM observation of Au NPs on r-TiO 2 (110) surfaces under O 2 at atmospheric pressure. Temperature series of HAADF-STEM images of TiO 2 -supported Au NPs during cooling under 10 5 Pa O 2 at (a) 400 °C, (b) 200 °C, (c) 100 °C and (d) RT. Unlike under Ar in Figure 1, particle ES evolves with temperature with the particle morphology becoming rounded at low temperature. (e)-(l) Temperature series of HAADF-STEM images over the temperature window where the octahedron-to-round morphology transition is observed. The temperature onset for this transition is situated at about 120 °C. (m) Quantitative follow-up of the extension of the {111}-and {100}-type segments delimiting the particle outline in projection as a function of temperature.

Figure 3 :

 3 Figure 3: In situ HAADF-STEM observation of Au NPs on r-TiO 2 (110) surfaces under H 2 at atmospheric pressure. Temperature series of HAADF-STEM images of TiO 2 -supported Au NPs during cooling under 10 5 Pa H 2 at (a) 400 °C, (b) 200 °C, (c) 100 °C and (d) RT. As under inert Ar gas in Figure 1, particle size and their truncated octahedron ES are stable down to RT.

Figure 3 2 nanorods acquired under 10 5

 325 Figure3shows a temperature series of HAADF-STEM images of Au NPs on r-TiO 2 nanorods acquired under 10 5 Pa H 2 just after the temperature has been successively set at 400 °C, 200 °C, 100 °C and RT. Similar to the in situ observation under Ar (Figure1), the NPs keep their initial sizes and morphologies as they are cooled to RT. The projected HAADF-STEM image of the Au NPs under H 2 are again compatible with truncated octahedron ES in dierent zone axes. The surface energy anisotropy ratio γ 100 /γ 111 , determined on four NPs, is equal to 1.07±0.04. This value is, within experimental errors, similar to the anisotropy ratio determined under Ar (1.10 ± 0.04).

2 and H 2 .

 22 γ int hkl allows us to construct the ESs (equilibrium shapes) of the supported Au NP under O All energetic values of the γ hkl and E ads , provided in the Supporting Information (Tables S1,S2 and S3), were calculated by performing accurate DFT calculations. Three low index Au (111), (100) and (110) surfaces were considered since they represent the main facets of stable truncated octahedron nanoparticles.

Figure 4 :

 4 Figure 4: Wul construction of Au NPs under oxygen environment at atmospheric pressure as a function of temperature from 200 °C to 17 °C.

Figure 5 :

 5 Figure 5: (a)The coverage of oxygen molecules on each facets of the Au NPs. (b)The revised interface energies of the Au NPs.(c) The facet fraction of the supported Au NPs.

Figure S4 :

 S4 Figure S4: Tested adsorption sites and adsorption congurations of oxygen molecule adsorbed on (a) Au(100), (b) Au(110), and (c) Au(111) surfaces.

Figure S5 .

 S5 Figure S5. Tested adsorption sites and adsorption congurations of H 2 molecule adsorbed on (a) Au(100), (b) Au(111), and (c) Au(110) surfaces.

Figure S6 .

 S6 Figure S6. Adsorption congurations of atomic oxygen on (a) Au(100), (b) Au(111) and (c) Au(110) surfaces

  However, environmental gas TEM observations of Au NPs under O 2 and H 2 have until now been limited to low pressure conditions of little relevance to catalytic reactions. The present works ambitions to bridge this pressure gap in the study of Au NP under O 2 and H 2 . Moreover, despite the successful use of MRS model in other systems, it is also the rst time that the equilibrium structure of Au NP under high O 2 and H 2 pressure is modeled and analyzed theoretically.

  Thus, for dissociative adsorption of O 2 by gold, no change in particle ES is expected during cooling to RT.Concerning H 2 adsorption on gold, DFT calculations show weak exothermic adsorption energies (TableS5). The predicted ES from the MSR model shows perfect truncated octahedron that does not evolve with pressure and temperature variations as observed experimentally. The excellent agreement between experiment and modelling is a strong indication that H 2 molecules do not adsorb readily on the Au NPs. The direct adsorption of H 2 on Au edge and corners atoms has been observed in Au NPs deposited on alumina.[START_REF] Bus | Hydrogen chemisorption on Al 2 O 3 -supported gold catalysts[END_REF] In this system, it has been shown that H 2 adsorption is limited by kinetic factors to23 to 100 °C with no spill-over to face sites. Similarly, DFT calculations have predicted that a single low-coordinated Au atom which can belong to defective extended surfaces or nanoparticles is sucient to dissociate H 2 . 31,71 The dissociation process involves a positive energy barrier but is spontaneous on Au atoms directly bonded to four other atoms. 71 In the present work, it is highly probable that when the NPs are cooled to RT the number of low-coordinated Au sites that can bind H 2 molecules is low. Consequently, the NPs initially retain their high-temperature Wul-Kaishev ES at all stages. In summary, by performing GCTEM on model supported Au NPs, we have provided the rst real space observations of O 2 and H 2 reactivity on nanosized gold at atmospheric pressure. By monitoring the change in particle ES upon gas exposure, we provide experimental evidence that Au NPs supported on rutile titania (110) surfaces can directly bind O 2 . Equilibrium shape reconstruction of Au NPs using the multiscale structure reconstruction modeling undertaken in conjunction with in situ observations of particle morphology brings fundamental atomic level insights into the absorption processes enabling the identication of the active sites at the particle surface for O 2 adsorption. Experiments and modelling show that Au(100) and (110) facets are active in O 2 adsorption but are inactive in H 2 adsorption.

Table S1 .

 S1 The DFT calculated surface energies of the Au(100), Au(110), and Au(111)

	under Vacuum.

Table S2 .

 S2 Adsorption energies (eV) of oxygen molecule on Au(100), (110), and

Acknowledgement J.N acknowledges funding from the French National Research Agency through the TOTEM project, Grant Number ANR-17-CE07-0031. The authors are also grateful to Région Ile-de-France for convention SESAME E1845 for the support of the JEOL JEM-ARM200F electron microscope installed at Paris Diderot University. J. M acknowledges the funding from the China Scholarship Council (File No. 201700260147) through the 2017-2019 Sino-French Cai Yuanpei Program.

Supporting Information Available

The following les are available free of charge. Graphical TOC Entry Some journals require a graphical entry for the Table of Contents. This should be laid out print ready so that the sizing of the text is correct.

Inside the tocentry environment, the font used is Helvetica 8 pt, as required by Journal of the American Chemical Society.

The surrounding frame is 9 cm by 3.5 cm, which is the maximum permitted for Journal of the American Chemical Society graphical table of content entries. The box will not resize if the content is too big: instead it will overow the edge of the box. This box and the associated title will always be printed on a separate page at the end of the document.