Maxime Guillotte
email: maxime.guillotte@univ-poitiers.fr

Julien Godet
email: julien.godet@univ-poitiers.fr

Laurent Pizzagalli
email: laurent.pizzagalli@univ-poitiers.fr

Molecular dynamics study of mechanical behavior of gold-silicon core-shell nanowires under cyclic loading

Keywords: reversible plasticity, shape restoration, Simulation, Amorphous silicon, Cyclic strain, composite materials, interface

Composite systems at nanoscale are promising to design new materials combining different interesting properties, like high strength and good plasticity. In the present work, we compare mechanical properties of gold-silicon core-shell and pure gold nanowires subjected to several tensile loadings and compressive unloadings at 300 K and 1 K, in order to evaluate the ability of a hard amorphous silicon shell to stabilize a crystalline gold core during cyclic loading. The hard shell appears involved in the quasi-reversible macroscopic plasticity of core-shell NW, by promoting twins, stacking fault and isolated full dislocations, due to a confinement effect of the core. In addition, our simulations revealed that the reversible behavior is progressively lost over cycles. The analysis of the amorphous shell shows that large strains combined with high strain rate prevent self-healing of the shell, leading to shell failure and to reversibility loss.

Introduction

New technologies often rely on multifunctional materials combining at least two interesting properties for applications, e.g. low density and high strength, high conductivity and high wearresistance, or high strength and good ductility [START_REF] Zhang | Approaching the ideal elastic strain limit in silicon nanowires[END_REF]. However, these properties are usually antagonistic or cannot co-exist. Using composite materials is one promising solution to design such multifunctional materials. It is particularly true at the nanoscale, where the ratios of surface on volume and/or interface on volume become predominant and allow to tune physical properties. For almost two decades, the fast development of experimental techniques for designing and testing nanostructures allows to reveal remarkable modifications of physical properties at nanoscale compared to their bulk counterparts, notably regarding mechanical aspects [START_REF] Uchic | Sample dimensions influence strength and crystal plasticity[END_REF]. The high strength sustained by such nanostructures has been thought to be used in composite materials to make them stronger [START_REF] Misra | Lengthscale-dependent deformation mechanisms in incoherent metallic multilayered composites[END_REF][START_REF] Zhu | Ultra-strength materials[END_REF][START_REF] Gu | Multiscale modeling of the elastic behavior of architectured and nanostructured cu-nb composite wires[END_REF].

Examples of interesting core-shell systems are Ti-TiN nanostructures, which exhibit high strength, high thermal stability, and good plasticity properties [START_REF] Zhang | Core-shell structured titanium-nitrogen alloys with high strength, high thermal stability and good plasticity[END_REF], core-shell nanowires (NWs) with increased Young modulus [START_REF] Chen | Determination of young's modulus of ultrathin nanomaterials[END_REF], or core-shell nanoparticles [START_REF] Zheng | Phase boundary effects on the mechanical deformation of core/shell cu/ag nanoparticles[END_REF][START_REF] Kilymis | Ductile deformation of core-shell si-sic nanoparticles controlled by shell thickness[END_REF]. Other more complex composite nanostructures, such as nanoporous gold covered by a thin layer of hard oxide, have been shown to limit coarsening and to enhance the materials strength [START_REF] Biener | ALD Functionalized Nanoporous Gold: Thermal Stability, Mechanical Properties, and Catalytic Activity[END_REF]. The latter effect has been recently investigated in coreshell Cu-Ni nanoporous metal by molecular dynamics simulation [START_REF] He | Deformation mechanisms and ductility enhancement in core-shell cu@ni nanoporous metals[END_REF]. The improvement of the mechanical properties seems partly linked to the presence of many interfaces, as in Cu-Nb composite nano-architectured materials [START_REF] Gu | Multiscale modeling of the elastic behavior of architectured and nanostructured cu-nb composite wires[END_REF]. Core-shell (Al -amorphous-Si) nanostructures have also been shown to be resistant to cyclic indentation [START_REF] Steck | Deformation and fatigue resistance of al/a-si core-shell nanostructures subjected to cyclic nanoindentation[END_REF] due to core confinement [START_REF] Fleming | The effects of confined core volume on the mechanical behavior of al/a-si coreshell nanostructures[END_REF]. Interestingly, restoration shape has also been reported experimentally in Ag -amorphous-SiO 2 NWs under bending [START_REF] Vlassov | Ants Lhmus, and Rnno Lhmus. Shape restoration effect in agsio2 coreshell nanowires[END_REF].

The goal of the present study is to better un-derstand the mechanical properties under cyclic loading of core-shell nanostructures combining high strength and good ductility. Core-shell NWs have been considered for their 1D shape, that allows an easier understanding of mechanical properties. Gold is used for the core for its good ductile property [START_REF] Ryu | Improved modified embeddedatom method potentials for gold and silicon[END_REF], while the shell is done with amorphous silicon taken as high strength material (at least at small dimensions [START_REF] Kang | Brittle and ductile fracture of semiconductor nanowires -molecular dynamics simulations[END_REF][START_REF] Östlund | Brittle-to-Ductile Transition in Uniaxial Compression of Silicon Pillars at Room Temperature[END_REF] and compared to gold).

In a previous work, we observed that such a hard amorphous shell is at the origin of plasticity confinement, resulting in an almost perfect elasto-plastic behaviour of core-shell NW, with shape conservation [START_REF] Godet | Uniform tensile elongation in au-si core-shell nanowires[END_REF]. Here, we examined whether plasticity confinement could be retained during cyclic loading, by using molecular dynamics (MD) simulations. We first considered a bare gold NW at 300 K as a reference. Then, we performed similar calculations on gold-silicon core-shell NWs at 300 K and 1 K to unravel the effects of both the amorphous shell and the temperature. The results are discussed in the last section of this paper.

Methods

Our work is based on MD simulations, carried out with the open-source code LAMMPS [START_REF] Plimpton | Fast parallel algorithms for short-range molecular dynamics[END_REF]. We used a Modified Embedded Atom Method (MEAM) potential that has been previously optimized to describe Au-Si interactions [START_REF] Godet | Uniform tensile elongation in au-si core-shell nanowires[END_REF]. The timestep used throughout the study is ∆t = 1 fs to ensure time integration stability. In order to model the Au-a-Si core-shell NW (AuSi NW), we first build an amorphous silicon cylinder, before replacing its core by crystalline gold. The amorphous silicon NW is created by randomly placing silicon atoms in a cylinder of principal axis ẑ, with a length L z = 203.5 Å, a radius of 120 Å and centered in x = 0 and y = 0. Periodic boundary conditions are used in ẑ direction. Using a conjugate gradient algorithm, we make a first rough energy minimization to move apart atoms which are initially too close from each other. A MD run is next performed at 2200 K in the NVT ensemble (constant number of atoms, volume and temperature) during 100 ps. The system is then quenched from 2200 K to 1 K during 250 ps and finally a conjugate gradient energy minimization is done. The radial distribution function of the amorphous silicon obtained by this method is consistent with previous studies [START_REF] France-Lanord | Atomistic amorphous/crystalline interface modelling for superlattices and core/shell nanowires[END_REF][START_REF] Khalid Laaziri | High resolution radial distribution function of pure amorphous silicon[END_REF][START_REF] Pedersen | Optimal atomic structure of amorphous silicon obtained from density functional theory calculations[END_REF] (Fig. 1). In a second step, a cylinder of principal axis ẑ, centered in x = 0 and y = 0, with a radius of ≈ 61 Å and a length L z along ẑ is defined. This volume is emptied and next filled with FCC gold atoms according to the following crystallographic orientations: x = [110], ŷ = [001], ẑ = [110]. The gold lattice parameter for this potential is a = 4.07 Å. A MD run is then performed at 300 K during 1 ns and the dimension along ẑ is allowed to change in order to decrease the residual stress σ zz along the ẑ axis to 0 GPa. This run was long enough to thermalize the system. The relaxed system is characterized by a length of 203.7 Å, a gold core radius R core of ≈ 60 Å and a silicon shell thickness e of ≈ 45 Å, which corresponds to 137101 Au atoms and 222256 Si atoms. In addition, we prepare a pure gold NW (Au-NW) containing as many atoms as in the core of the AuSi-NW, with the same crystallographic orientations. After ẑ-relaxation at 300 K during 1 ns, final length and radius R Au were 205 Å and ≈ 59 Å respectively. This run was long enough to thermalize the system. We performed NVT tension-compression deformation cycles on AuSi-NW at 300 K and 1 K as well as on Au-NW at 300 K, by tuning the box dimension along ẑ. The engineering strain rate is equal to 10 8 s -1 , which is typical of MD simulations. The maximum deformation reached in our simulations corresponds to an engineering strain of 30 %. To determine the engineering stress, the computed load is divided by the wires section before deformation (π(R core + e) 2 and π(R Au ) 2 for the AuSi-NW and the Au-NW, respectively).

Results

Au-NW at 300K

We first focus on the Au-NW, since its behavior during cyclic deformation will serve as a reference for the AuSi-NW study. Fig. 2 shows the engineering stress versus the engineering strain for four tension-compression cycles on Au-NW at 300 K. During the first tensile loading, the Au-NW is elastically deformed until = 5.0 % with a Young modulus of 74 GPa. The first plastic event is an avalanche of spatially and temporally correlated dislocations [START_REF] Godet | Uniform tensile elongation in au-si core-shell nanowires[END_REF] leading to the nucleation and growth of a twin as already shown in [START_REF] Liang | Atomistic simulations reveal shape memory of fcc metal nanowires[END_REF][START_REF] Li | Surface rebound of relativistic dislocations directly and efficiently initiates deformation twinning[END_REF], and the nucleation of a lonely dislocation leading to an isolated stacking fault. We carried out a dislocation analysis with the visualization and analysis software Ovito [START_REF] Stukowski | Visualization and analysis of atomistic simulation data with ovito -the open visualization tool[END_REF] for a range of deformations corresponding to the first drop of the black curve in Fig. 2. Three Shockley partial dislocations (propagating in the [112] and [112] directions in three consecutive (111) planes) are nucleated one after another leading to twin formation. A snapshot of the microstructure just after the first stress drop at = 5.2 % is shown in Fig. 3(b). Up to = 14.6 %, the plastic deformation mainly corresponds to the growing of this twin (Fig. 3 (b-c)), by the propagation of leading partial dislocations on the upper twin boundary plane and along [112] (Fig. 3(bc)). Between = 14.6 % and = 30 %, the plastic deformation is next associated with full slip (propagation of dissociated perfect dislocations) on the upper twin boundary plane (Fig. 3(c-d)).

Experimentally, in thin films, the transition from twinning to full dislocation slips has been attributed to the increase of the films thickness [START_REF] Oh | In situ tem straining of single crystal au films on polyimide: Change of deformation mechanisms at the nanoscale[END_REF]. In NWs, however, twin formation is not necessarily followed by dislocations slip, but rather by twinning on conjugate slip planes [START_REF] Lee | Reversible cyclic deformation mechanism of gold nanowires by twinning-detwinning transition evidenced from in situ TEM[END_REF]. In our calculations, we observed a transition from twinning to full slip, in agreement with previous atomistic simulations [START_REF] Liang | Atomistic simulations reveal shape memory of fcc metal nanowires[END_REF][START_REF] Rezaei | Pseudoelasticity and shape memory effects in cylindrical FCC metal nanowires[END_REF]. This difference between experiments and simulations might be partly explained by the crystal re-orientation caused by twinning and the use of periodic boundary conditions in simulations, whereas experimentally NW ends are clamped. In the simulation, as the twin size increases, the stress axis slightly tilts, thus modifying the Schmid factors on leading (from 0.471 to 0.498) and on trailing partial dislocations (from 0.236 to 0.242) before and after the twin formation. In addition, the growth of a single twin leads to ledges formation at the intersection of surface and twin boundaries, which act as local stress concentrators for dislocations nucleation [START_REF] Brochard | Elastic limit for surface step dislocation nucleation in facecentered cubic metals: Temperature and step height dependence[END_REF]. These two effects enhance the nucleation probability of leading, but also trailing partial dislocations. Once a trailing dislocation is nucleated, a sharp step is created on the surface, which increases the stress concentration on a unique slip plane leading to the emission of several dissociated perfect dislocations in a single plane. This mechanism of twinning followed by full slip during tensile loading, greatly affects the morphology of the Au-NW, as shown in Fig. 3(d).

During unloading, the Au-NW is elastically deformed between = 30 % and = 24.4 %. From 24.4 % to 14.2 %, detwinning occurs by the nucleation and migration of partial dislocations in the lower twin boundary plane (Fig. 3(d)). After detwinning, the plastic deformation is associated to multiple dislocations slip in the same plane between = 14.2 % and = 0 % (Fig. 3(e-f)). The sequence of twinning followed by full slip in traction, and detwinning followed by full slip in compression, is related to different Schmid factors on the leading and trailing partial dislocations for traction and compression, as explained in previous works [START_REF] Park | Deformation of fcc nanowires by twinning and slip[END_REF][START_REF] Seo | Superplastic deformation of defect-free au nanowires via coherent twin propagation[END_REF][START_REF] Christopher | Atomistic simulations and continuum modeling of dislocation nucleation and strength in gold nanowires[END_REF]. Note that twinning and dislocation slips can also depend on the surface roughness of NWs as already shown in Ref. [START_REF] Park | Deformation of fcc nanowires by twinning and slip[END_REF].

In compression, serrated stress variations have a larger amplitude than in tension, due to the cas- cading nucleation of full dislocations. Conversely, under tensile loading, dislocations tend to be emitted one by one (Fig. 2). Moreover, full dislocations migrate in the same slip plane, and each nucleation event increases the surface step height under tensile loading (Fig. 3(e-f)). Since stress concentration increases with step height, lower applied stress are required for nucleating additional dislocations. Instead, in compression, each nucleated dislocation decreases the step height (Fig. 3(e-f)), and nucleating further dislocations requires larger applied stresses.

After full unloading, the wire globally recovers its shape and all the crystallographic defaults in the bulk (stacking faults and twins) are healed (Fig. 3(f)). However, the NW surface is now rougher and the wire at = 0 % is in a compressive state (σ zz ≈ -2 GPa). This can be explained by the fact that there were less dislocation slip events during unloading than loading stage, and thus the NW stretching is not fully compensated over the four cycles.

During the second cycle, the Au-NW is elastically deformed from = 0 % to = 5.3 % and the new mechanical equilibrium (σ zz = 0) is reached at = 2.4 %. For the rest of this cycle and during the subsequent cycles, twinning is followed by full dislocation slips during tensile loading, while during compressive loading, detwinning is followed by dislocation slips. For second and subsequent cycles, yield points are quasi equivalent and the stressstrain curves are almost superimposed. Note that the yield stresses are slightly lower than the value obtained for the first cycle, essentially because surface defects created during the first cycle facilitate the plastic deformation for the next cycles compared to the initial pristine system. We also observed that the surface roughness increases with cycles, at = 0 % (Fig. 3(f,h,j,l)) and at = 30 % (Fig. 3(d,g,i,k)).

AuSi-NW at 300K

Fig. 4 shows engineering stress versus engineering strain curves for six tension-compression cycles and the core-shell AuSi-NW, at 300 K. During the first tensile loading, the AuSi-NW is elastically deformed until = 5.8 %, with a Young modulus equal to 90 GPa, larger than the Au-NW one because of the silicon shell. The stress-strain curves do not show the sawtooth behavior (Fig. 4) as obtained for Au-NW (Fig. 2). Indeed, from 5.8 % to 10.1 %, only partial dislocations are nucleated, with no dislocations avalanche, thus limiting stress drops. Also, there is no spatial correlation between two successive dislocations, as well as no twin formation (Fig. 5(b). The silicon shell is the main factor limiting dislocation cascades, and delaying twin growth. At strains greater than = 10.1 %, twins begin to grow, accompanied by the nucleation and propagation of isolated full and partial dislocations (Fig. 5(b-c)). Due to the silicon shell confinement, dislocations are homogeneously distributed all along the interface and the NW core morphology is preserved. However, roughness of the shell surface increases compared to the initial structure. During the first unloading, the AuSi-NW is elastically deformed between = 30 % and = 25 %. From = 25 % to = 0 %, similar plastic deformation mechanisms than during loading are observed, i.e. nucleation of partial dislocations in adjacent planes leading to detwinning, and nucleation and propagation of isolated partial and full dislocations (Fig. 5(c-d)). At the end of the first unloading, all bulk defects in the core have been healed. However, both the shell surface and the core-shell inter-face appear slightly rougher compared to the initial state. Furthermore, a limited interdiffusion is occurring at the interface, in agreement with a previous work on Au-Ni core-shell NW [START_REF] Abdolrahim | The mechanical response of core-shell structures for nanoporous metallic materials[END_REF]. Note that the AuSi-NW is in a compressive state at = 0 % (σ zz ≈ -3.7 GPa), like Au-NW after the first cycle, for similar reasons.

During the second cycle, the NW is elastically deformed until = 6.7 % and reach mechanical equilibrium at = 4.9 % (Fig. 4). The stress-strain curve is different from that of the first cycle, with a lower stress maximum occurring at a larger strain. The flow stress is also slightly lower although plasticity mechanisms are similar. However, the unloading part closely resembles the one in the first cycle, both in terms of stress variations and deformation mechanisms. Stress-strain curves are almost superimposed for the following cycles, except regarding the stress reached at maximum strain. In fact, the latter clearly decreases as the numbers of cycles grow. Examining in details the morphology of the AuSi NW during cycles, one can observe the formation of trenches at the shell surface (Fig. 5(c,e,g,i)), deepening as a function of cycles. Those are only partially healed during the unloading stages (Fig. 5(d,f,h,j)). Furthermore, cyclic deformation tends to slowly change the core morphology, in an apparent relation with shell trenches.

We also tested the influence of the shell thickness on the mechanical properties of the AuSi-NW. We found that the stress-strain curve and the deformation mechanisms obtained with a shell as thin as 1 nm, are similar to those obtained with the larger shell of 4.5 nm (See Supplementary Material). This suggests that the results presented in this work are not specific to systems with thick shells.

AuSi-NW at 1K

In order to study the influence of the temperature, four tension-compression cycles have been performed for the AuSi-NW, at 1 K. The calculated stress-strain curves are shown in Fig. 6. Compared to the 300 K test, the elastic deformation stage lasts until = 5.8 %, with a Young modulus of 110 GPa. This is higher than at 300 K, in agreement with the well accepted fact that the materials stiffness is greater at 0 K than at 300 K. For the four tensile loadings, identified plastic events include a few isolated partial dislocations and a few isolated dissociated perfect dislocations, and mainly the growth of a unique twin (Fig. 7(a,c,e,g)). At the highest strains, this large twin induces a necking of the gold core, correlated to the formation of deep and large trenches in the silicon shell. However, the core necking is lighter than the one occurring at the upper twin boundary in the bare Au-NW.

During unloading stages, the activation of plastic deformation mechanisms like detwinning allows for recovering a gold core containing no or few residual defects (Fig. 7(b,d,f,h)). The opened trenches in the shell are also partially bridged. However, both the shell and the core morphology keep stigma associated to the cyclic deformation.

Discussion

The singularity of the first loading

For all tested systems, we find that the stress variations monitored during the first loading are different than those obtained for the following cycles. In particular, larger maximum stress values are reached (See Fig. 2, 4, and6). This is especially true for core-shell systems. We tentatively attribute this result to structural changes between the initial state and the beginning of the second loading. For instance, the starting configuration for the pure gold NW is alike a perfect cylinder without any surfaces defects. After the first cycle, the surface is more rough, with the presence of atomic steps left by dislocations. It is well known that such surface defects are efficient stress concentrators [START_REF] Brochard | Elastic limit for surface step dislocation nucleation in facecentered cubic metals: Temperature and step height dependence[END_REF], favoring the nucleation of dislocations at a lower applied stress. For the core-shell NW, a low density of defects is present at the Au-Si interface for the initial structure, due to long annealing and relaxation. After the first cycle, the amount of interface defects has substantially increased, as well as the interface roughness (Compare Fig. 5-a and Fig. 5-d). Again, it facilitates the nucleation of the first dislocations, which then happened at a lower applied stress than during the first cycle.

The singularity of the first loading compared to the followings ones is then related to the use of initial configurations with no or few structural defects. It is also probably an indirect consequence of the high strain rates used in molecular dynamics simulations. In fact, the short simulation time does not allow for an efficient thermal annealing of the structural defects left by plastic deformation. The influence of thermal annealing can be easily estimated by comparing the stress -strain curves and configurations at 300 K (Figs. 45) and 1 K (Figs. 67), respectively. At 300 K, the maximum stress drops by about 1.1 GPa between the first and second cycle, whereas is is about 2.1 GPa at 1 K. The interface roughness is also more pronounced at 1 K than at 300 K. These observations suggest that in experiments, with characteristic times of the order of the minute or the hour, the difference between the first cycle and the following ones should be greatly reduced compared to the present work.

Reversibility of the Au NW plastic deformation

Putting aside the first cycle and focusing only on the following ones, the stress -strain curves for the Au NW seem to indicate that the plastic deformation is reversible (Fig. 2). This feature is related to the fact that twinning is the favored mode of deformation at low strains. This mechanism is at the origin of the well documented pseudoelastic behavior of metallic nanowires with square and rhombic cross sections [START_REF] Liang | Atomistic simulations reveal shape memory of fcc metal nanowires[END_REF][START_REF] Lee | Reversible cyclic deformation mechanism of gold nanowires by twinning-detwinning transition evidenced from in situ TEM[END_REF][START_REF] Park | Shape memory and pseudoelasticity in metal nanowires[END_REF]. However, in the present work, twinning is always accompanied with full dislocation slip events, probably because of the size and circular cross-section of the NW. Dislocation slip leads to a significant shearing of the NW which is only partially healed during unloading, conversely to twinning. This is confirmed by the seemingly roughness increase of the Au NW surface during cycles (Fig. 4). The cyclic deformation of the Au NW is then better described as quasi-reversible. Note that only a limited number of cycles have been performed in this work, and it is not clear whether additional cycles will continue to change the NW structure, or if a stable state will be reached.

Reversibility of the AuSi NW plastic deformation

We now discuss the reversibility of the AuSi NW plastic deformation under cyclic deformation. At both 300 K (Fig. 4) and 1 K (Fig. 6), the stress -strain curves show a significant overlap for the second and next cycles. However the reversible character is less prominent than for bare Au NW, in particular the stress values at high strain now decrease as a function of cycles. For further insights, we plot in Fig. 8 the axial stress from the NW center to its surface, at = 30 %, 300 K and for each cycle. The curves are obtained by computing the mean atomic stress concentric hollow cylinders with equal thickness at different radii. When integrated over the whole NW (dashed curves in Fig. 8), the stress values at 30 % as reported in Fig. 4 are recovered.

Examining first the core part (for distance lower than 60 Å), there is no clear trend regarding the stress variations between each cycle. Those latter are due to different core structures after each loading. Conversely, in the shell part, the axial stress decreases as a function of cycles. It is noteworthy that the more the radial distance increases, the more the weight of the stress contribution to the mean stress is important, since the number of atoms contained in a cylinder is proportional to its squared radius. This is evidenced by integrating the axial stress on a cylinder of increasing radius (Fig. 8, dashed curves). Therefore, the stress decrease as a function of cycles observed in Figs. 4,6 can be mainly attributed to the shell.

The analysis of the evolution of the amorphous shell structure over cycles reveals the progressive formation of trenches on the shell surface, that weakens the shell mechanical resistance to loading, and is certainly the cause of the aforementioned stress decrease. The consequences are more dramatic at 1 K, with a localized breaking of the shell (Fig. 7-a,c,e,g and 9-a,b). These figures show a spatial correlation between the large twin formed in the core and those trenches, at = 30 %. During tensile loading, the formation of the twin leads to crystal rotation and core necking, generating shear strain in the shell. In fact, examining the atomic shear strain just before the rupture of the shell suggests that the largest values correspond to twin regions in the vicinity of the trenches (Fig. 9-c). The combination of this shear strain with the applied tensile stress induces a mixed loading on the shell under mode 1 and 2 leading to the expansion of the trenches, and potentially to shell fracture. This analysis is in agreement with the brittle nature of silicon, in particular at low temperature [START_REF] Korte | Deformation of silicon -insights from microcompression testing at 25 -500 degree c[END_REF][START_REF] Rabier | Plastic deformation of silicon between 20c and 425c[END_REF].

All these elements indicate that the deformation of AuSi-NW under cyclic loading is not perfectly reversible, and that the irreversible part can be essentially ascribed to the shell. We propose the following scenario. During loading, there is a stress build-up on silicon atoms which is only partially relaxed by small atomic displacements, typical of amorphous plastic deformation. At high strains, a more efficient stress relaxation is obtained by the breaking of silicon bonds, leading to the formation of trenches. The shell structure integrity is not fully recovered during unloading (Fig. f,h), resulting in the progressive opening of the trenches.

Influence of various parameters on shell weakening

In this part we examine the effect of different parameters on the weakening of the Si shell, and thus on the reversibility of the AuSi plastic deformation. Two different temperatures, 300 K and 1 K, have been tested in this work. Comparing stress -strain curves, in Figs. 4 and6, indicates that the reversible character is stronger at 300 K than at 1 K. In fact, the stress at 30% decrease by 0.5 GPa from the first to the fourth cycle in the former case, whereas the stress reduction is about 1.1 GPa in the latter case. This suggests a thermally activated recovery of the Si shell during the consecutive extension-compression cycles. This conclusion is supported by the examination of the shell structure and morphology. Significant cracks already occurs during the 2 nd cycle at 1 K (Fig. 7), although at 300 K the shell is not fully broken even after the 6 th cycle (Fig. 5). Overall, the morphology of the shell is more preserved at 300 K than at 1 K.

Another important factor is the strain rate inherently used in MD simulations, which is typically 10 11 times higher than in experiments. Such a large difference makes it almost impossible to use the same procedure to study the effect. In fact, decreasing the strain rate to 10 7 s -1 or even lower would require a large increase of computational resources, with a likely unnoticeable influence. To circumvent this issue, we carried out additional NVT calculations at fixed strain (30%), which then mimic a deformation test with an infinitely slow strain rate. The Fig. 10 shows the results obtained at 300 K over a MD time of 3 ns, for two structures obtained by cyclic deformation at a strain of 30% and at cycles one and six. For the first (sixth) cycle, the axial stress decreases from 2.9 GPa to 2 GPa (from 2.2 GPa to 1.5 GPa). This result suggests that significantly reduced stresses would be obtained if low strain rates could be used in MD simulations. At 10 8 s -1 , there is not enough time for thermally acti-vated stress relaxation mechanisms to occur, which leads to stress build-up in localized regions, and finally trenches formation and cracks propagation. The stress variations in Fig. 10 have been fitted with exponentially decreasing functions (fitting parameters reported in the figure), to estimate the residual axial stress at times greater than 3 ns. Choosing 1 s as a time typical for experimental deformation of nanowires, the axial stress is extrapolated to 1.63 GPa for the first, and 1.38 GPa for the sixth cycle. The difference between these two values is only 0.25 GPa, to be compared to the initial stress difference of 0.75 GPa between the first and the sixth cycles (at 30% in Fig. 4). Therefore, the extrapolation to times comparable to experiments suggests a significant improvement of the amorphous shell and by consequence of the reversibility of the plastic deformation of the coreshell nanowire during cyclic deformation.

At last, the influence of the maximum strain attained during tensile loading is discussed. Cyclic deformation of the AuSi NW between 0% and 9% have been carried out at 1 K, for three cycles. This specific strain value has been selected since it corresponds to the strain threshold above which twinning is occurring. The stress -strain curves and snapshots of 9% strained configurations are shown in Supplementary Materials. They reveal a greatly improved reversibility compared to the case with 30% maximum strain. In fact, there is now no stress decrease at maximum strain after the first cycle. Also the shell appears to be completely preserved during the cyclic deformation. Consequently, it can continue to confine and spatially delocalize the plastic deformation in the gold core with maximal efficiency.

Conclusion

In this work, we focused on the mechanical properties of gold NWs and Au-Si core-shell NWs during cyclic loading, in order to investigate the influence of a hard amorphous silicon shell on gold NW. We first observed an increase of the yield stress due to the presence of the hard shell, compared to pure gold NW. We also found that a hard shell can confine the plasticity in the core, leading to shape conservation at deformation as large as 30% in tension over several cycles. The presence of the shell favors twins, isolated partials and isolated full dislocations over localized full dislocations. This allows for a seemingly macroscopic reversible plas-ticity, with the presence of interface defects. Usually, the mechanisms of reversible plasticity is based on twinning and detwinning that allows for an homogeneous distribution of plasticity all along the NW length. Here, the hard amorphous shell favors not only twinning but also an homogeneous distribution of full dislocations. Therefore, a larger amount of plastic strain could be stored in these composite core-shell NWs than in Au NWs, which should lead to reversible plasticity at larger strain. However, we also found that the reversible plasticity effect is progressively lost at 300 K over cycles and completely lost at 1 K after the first loading. The analysis of thermal annealing of strained core-shell NW at 300 K on a long timescale, suggest that the high strain rate in MD simulations induces a limited self-healing of the amorphous shell, thus leading to defects storage during cyclic loading. These defects appear as precursors of cracks and shell rupture. Additional analyses also show that shell degradation can be greatly reduced if the maximum loading strain remains below the threshold above which twinning is occurring. Our investigations imply that optimizing these different factors (temperature, strain rate, maximum loading strain) could enhance the reversibility of the plastic deformation during cyclic tests. This obviously suggests potential directions for future researches.

Figure 1 :

 1 Figure 1: Radial distribution function of an experimental amorphous silicon [21] (red curve) ; Radial distribution function of our simulated amorphous silicon after the thermalization at 300 K (black curve).

Figure 2 :

 2 Figure 2: Engineering stress -engineering strain curves for 4 tension-compression cycles done at 300 K for Au-NW.

Figure 3 :

 3 Figure 3: Evolution of Au-NW structure (cross-section views) during the first cycle (a-f). Structures at the end of the second, third and fourth tensile loading (g,i,k respectively) and structures at the end of the second, third and fourth compressive loading (h,j,l respectively). Purple, yellow, red and grey spheres show respectively Au atoms of the parent crystal (FCC environment), Au atoms of the twinned crystal (FCC environment), Au atoms in HCP environment and Au atoms in unknown environment. The full (dashed) line arrows indicate the direction of twin boundaries (dislocation) motion.

Figure 4 :

 4 Figure 4: AuSi-NW engineering stress -engineering strain curves for 6 tension-compression cycles done at 300 K.

Figure 5 :

 5 Figure 5: Evolution of AuSi-NW structure (cross-section views) during the first cycle (a-d). Structures at the end of the second, fourth and sixth tensile loading (e,g,i respectively) and structures at the end of the second, fourth and sixth compressive loading (f,h,j respectively). Purple, yellow, red, grey and green spheres are respectively Au atoms of the parent crystal (FCC environment), Au atoms of the twinned crystal (FCC environment), Au atoms of HCP environment, Au atoms of unknown environment and silicon atoms. The red arrows indicate the location of trenches which form on shell surface during cycles.

Figure 6 :

 6 Figure 6: AuSi-NW engineering stress -engineering strain curves of 4 cycles tension-compression performed at 1 K.

Figure 7 :

 7 Figure 7: At 1 K: Microstructures at the end of the first, second, third and fourth tensile loading (a,c,e and g respectively) and microstructures at the end of the first, second, third and fourth compressive loading (b,d,f and h respectively) (See Fig. 5 for color coding).

Figure 8 :

 8 Figure 8: Average axial stress of the AuSi-NW at = 30 % during the cyclic deformation at 300 K as a function of the radial distance for each cycle (full lines). Stress integrated on a cylinder of radius r and divided by the area of the initial cross section S 0 ( r 0 2πrσdr/S 0 ) (dashed lines).

Figure 9 :

 9 Figure 9: (a-b) 3D representations of the AuSi-NW at the end of the fourth tensile loading at 1 K, (a: full Au-Si-NW, b: gold core only). A cross-sectional view is shown in Fig. 7(g). The color coding is the same than in Fig. 3, 5 and 7. (c) Atomic shear strain associated to core gold atoms in AuSi-NW for the same system.

Figure 10 :

 10 Figure 10: Evolution of the AuSi Nw engineering stress during thermal annealing, for two initial configurations corresponding to NW pulled at 30%, at the end of the first and the sixth tensile loading (black and red curves respectively), with their fitting curves (green and blue curves respectively).

Acknowledgements

This work pertains to the French Government program Investissements d'Avenir (LABEX IN-TERACTIFS, reference ANR-11-LABX-0017-01). Computations have been performed on the supercomputer facilities of the Mésocentre de calcul Poitou-Charentes.


Mohammad Javad Shojaei

Antonio Rodríguez De Castro

Yves Méheust

Dr Nima Shokri
email: nima.shokri@manchester.ac.uk

Dynamics of foam flow in a rock fracture: Effects of aperture variation on apparent shear viscosity and bubble morphology

Keywords: Fractured media, Foam flow, Aperture heterogeneity, Shear rheology, Bubble morphology

There has recently been renewed interest in understanding the physics of foam flow in permeable media. As for Newtonian flows in fractures, the heterogeneity of local apertures in natural fractures is expected to strongly impact the spatial distribution of foam flow.

Although several experimental studies have been previously performed to study foam flow in fractured media, none of them has specifically addressed that impact for parallel flow in a realistic fracture geometry and its consequences for the foam's in situ shear viscosity and bubble morphologies. To do so, a comprehensive series of single-phase experiments have been performed by injecting pre-generated foams with six different qualities at a constant flow rate through a replica of a Vosges sandstone fracture of well-characterized aperture map.

These measurements were compared to measurements obtained in a Hele-Shaw (i.e., smooth) fracture of identical hydraulic aperture. The results show that fracture wall roughness strongly increases the foam's apparent viscosity and shear rate. Moreover, foam bubbles traveling in regions of larger aperture exhibit larger velocity, size, a higher coarsening rate, and are subjected to a higher shear rate. This study also presents the first in situ measurement of foam bubbles velocities in fracture geometry, and provides hints towards measuring the in situ rheology of foam in a rough fracture from the velocity maps, for various imposed mean flow rates. These findings echo the necessity of considering fracture wall when predicting the pressure drop through the fracture and the effective viscosity, as well as in situ rheology, of the foam.

Introduction

Fracture media are present in a wide range of geological media used for industrial applications including CO 2 sequestration, subsurface/soil remediation, and Enhanced Oil Recovery (EOR) [START_REF] Wu | Multiphase Fluid Flow in Porous and Fractured Reservoirs[END_REF][START_REF] Chen | Pore-scale simulation of density-driven convection in fractured porous media during geological CO2 sequestration[END_REF][START_REF] Cicek | Compositional and non-isothermal simulation of CO2 sequestration in naturally fractured reservoirs/coalbeds: Development and verification of the model[END_REF][START_REF] Chen | Experimental demonstration of contaminant removal from fractured rock by boiling[END_REF][START_REF] Singhal | Applied hydrogeology of fractured rocks[END_REF]. In these applications, when injecting a displacing fluid in a reservoir containing fractures, channeling and preferential flow paths occur, which leads to low sweep efficiency [START_REF] Tsang | Flow channeling in heterogeneous fractured rocks[END_REF][START_REF] De Wit | Viscous fingering in periodically heterogeneous porous media. II. Numerical simulations[END_REF][START_REF] Al-Wahaibi | First-contact-miscible and multicontact-miscible gas injection within a channeling heterogeneity system[END_REF]. The use of foams, which consist of gas bubbles separated by thin liquid films denoted as lamellae [START_REF] Osei-Bonsu | Foam stability in the presence and absence of hydrocarbons: From bubble-to bulk-scale[END_REF][START_REF] Hirasaki | Mechanisms of foam flow in porous media: apparent viscosity in smooth capillaries[END_REF][START_REF] Shojaei | Foam stability influenced by displaced fluids and by pore size of porous media[END_REF], has proven to be a potential solution to overcome these issues [START_REF] Abbaszadeh | Single-Well Simulation Study of Foam EOR in Gas-Cap Oil of the Naturally-Fractured Cantarell Field[END_REF][START_REF] Kvaestad | CO2-foaming agent retention in fractured chalk models: Experiments and simulations[END_REF][START_REF] Johansen | An Experimental Study of Foam Flow in Fractured Systems of Increasing Size[END_REF][START_REF] Telmadarreie | Evaluating the Potential of CO2 Foam and CO2 Polymer Enhanced Foam for Heavy Oil Recovery in Fractured Reservoirs: Pore-Scale and Core-Scale Studies[END_REF]. In particular, a high viscous pressure drop is produced through the foamsaturated fracture, diverting the flow to less permeable regions and thus improving oil recovery [START_REF] Haugen | Experimental study of foam flow in fractured oilwet limestone for enhanced oil recovery[END_REF][START_REF] Farajzadeh | Foam assisted gas-oil gravity drainage in naturallyfractured reservoirs[END_REF][START_REF] Fernø | Mobility control during CO 2 EOR in fractured carbonates using foam: Laboratory evaluation and numerical simulations[END_REF][START_REF] Zuta | Transport of CO2-foaming agents during CO2-foam processes in fractured chalk rock[END_REF].

Foams are generated by three main mechanisms: snap-off, lamellae division, and leavebehind [START_REF] Ransohoff | Mechanisms of foam generation in glass-bead packs[END_REF][START_REF] Shojaei | Foam Flow Investigation in 3D-Printed Porous Media: Fingering and Gravitational Effects[END_REF], and can exist in either bulk or confined form [START_REF] Andrianov | Immiscible foam for enhancing oil recovery: bulk and porous media experiments[END_REF]. "Bulk foams" are foams with individual bubbles which are considerably smaller than the characteristic length scale of the porous medium [START_REF] Aronson | The influence of disjoining pressure on foam stability and flow in porous media[END_REF]. This type of foams can be classified into ball (or wet) foam and polyhedral (or dry) foam, depending on the foam bubble shape [START_REF] Kovscek | Fundamentals of foam transport in porous media[END_REF]. Polyhedral foam bubbles are more common at high gas volume fractions, while spherical foam bubbles form at low gas volume fraction [START_REF] Kovscek | Fundamentals of foam transport in porous media[END_REF][START_REF] Verbist | The foam drainage equation[END_REF]. On the other hand, confined foams are produced when the characteristic length scale of individual bubbles is of the same order of magnitude or greater than, the characteristic length scale of the porous medium [START_REF] Kovscek | Engineering, Foam flow through a transparent rough-walled rock fracture[END_REF].

Predicting and controlling foam coalescence is a key factor for successful foam flooding projects. Chambers and Radke [START_REF] Chambers | Interfacial phenomena in petroleum recovery[END_REF] proposed two main mechanisms of foam coalescence.

The first one is capillary-suction, which results from the Lamellae's expansion and contraction to conserve mass as they go out through constrictions in the porous medium. The lamellae thin and thicken in a sequence of squeezing-stretching and draining-filling events. This oscillation of the lamella thickness has a wider oscillation range at higher gas flow rate and larger pore-body to pore-throat aspect ratio [START_REF] Jiménez | Dynamic stability of foam lamellae flowing through a periodically constricted pore[END_REF], and can lead to lamella rupture. In this respect, Jimenez and Radke [START_REF] Jiménez | Dynamic stability of foam lamellae flowing through a periodically constricted pore[END_REF] stated that the capillary suction pressure increases with porebody to pore-throat aspect ratio. The other mechanism stems from Ostwald ripening and occurs due to differences in gas pressure between neighboring bubbles of different sizes, due to Laplace law, which stipulates that the pressure drop across the interface is all the larger as the curvature ratio of the interface is smaller. That pressure difference results in gas diffusion from the smaller to the larger bubbles, and progressive vanishing of the smaller bubbles.

According to Von Neumann's law [START_REF] Neumann | Metal interfaces[END_REF], the rate of change in bubble surface area depends on the number of neighbor bubbles, the surface tension of the liquid and the permeability of the porous medium. It is worth mentioning that whereas foam coarsening was extensively investigated in bulk systems, there is less evidence of this phenomenon in confined geometry [START_REF] Jones | Foam flow in a model porous medium: I. The effect of foam coarsening[END_REF], probably because the involved diffusive time scales are larger than the characteristic time of most laboratory experiments.

Foams can serve as efficient displacing fluids in porous media due to their specific rheological properties. A concept frequently used to characterize the rheology of foams is their apparent (or effective) viscosity, (Pa.s). Knowing the foam's apparent injection 𝜇 app velocity (or specific discharge) q (m/s), defined as the ratio of the total volumetric flow rate of the foam to the cross-sectional area of the fracture , with w being the width of 𝑆 = 𝑤 × ℎ m the fracture and (m) its mechanic (that is, arithmetic mean) aperture, the pressure drop ℎ m ∆𝑃 through the fracture of length L (m), and the permeability K (m 2 ) of fracture, the (𝑃𝑎) apparent viscosity is the viscosity that must be used in order to be able to write a Darcy's law for the foam flow [START_REF] Neumann | Metal interfaces[END_REF]:

(1)

𝜇 app = 𝐾∆𝑃 𝑞𝐿 = ℎ 2 ∆𝑃 12𝑞𝐿
In the preceding equation, (m) is the constant aperture of a smooth fracture ℎ = (12𝐾) 1/2 (parallel plate) generating the same pressure drop as the considered real fracture when a Newtonian fluid is injected at the same flow rate.

It is important to note that is not an intrinsic property of the foam: its value is usually 𝜇 app observed to decrease as q increases, which indicates a shear-thinning behaviour of the foam. Indeed, the foam bubbles tend to agglomerate at low flow rates, which impedes their relative movement, but are aligned and deformed in the direction of flow when the flow rate increases, resulting in a decrease in viscosity. Note however that the in-situ rheology of a foam in a porous medium whose characteristic pore size is on the same order as the typical foam bubble size, is not necessarily the same as that of the bulk flow of the same foam, since in the former case most friction occurs between the liquid films and the solid walls, while in the latter case viscous dissipation occurs mostly within liquid films between bubbles. In this regard, Hirasaki and Lawson [START_REF] Hirasaki | Mechanisms of foam flow in porous media: apparent viscosity in smooth capillaries[END_REF] showed that the apparent viscosity of a foam is the sum of three contributions: (1) the flow resistance of slugs of liquid between gas bubbles (which results from viscous friction at the solid walls), (2) the flow resistance due to surface tension gradient between the front and rear of foam bubbles (3) and the resistance to deformation of foam bubbles as they pass through the pores. They considered both bulk foam flows and the flow of foam in porous media. They also showed that foam texture (i.e. the number of bubbles per unit area) is a key parameter controlling . Some studies showed that foam 𝜇 app behaves as a yield stress fluid [START_REF] Rouyer | Is the yield stress of aqueous foam a well-defined quantity?[END_REF][START_REF] Lindner | Viscous fingering in a yield stress fluid[END_REF][START_REF] Gardiner | Yield stress measurements of aqueous foams in the dry limit[END_REF]. After the stress passes a threshold, foam bubbles rearrange and flow like a viscous non-Newtonian fluid. Princen [START_REF] Princen | Rheology of foams and highly concentrated emulsions: I. Elastic properties and yield stress of a cylindrical model system[END_REF] and Hohler and Cohen-Addad [START_REF] Höhler | Rheology of liquid foam[END_REF] showed that the yield stress is inversely proportional to the bubble size in bulk and bubble scale.

The effects of aperture variations on the creeping (that is, laminar with ) flow through 𝑅𝑒 ≪ 1 a rough fracture, in particular its departure from the ideal model of the flow through two perfectly smooth parallel palates [START_REF] Brown | Fluid flow through rock joints: the effect of surface roughness[END_REF], have been addressed extensively in the literature of the last 35 years [START_REF] Brown | Fluid flow through rock joints: the effect of surface roughness[END_REF][START_REF] Méheust | Flow enhancement of a rough fracture[END_REF][START_REF] Méheust | Geometrical heterogeneities and permeability anisotropy of rough fractures[END_REF]. Spatial variations of aperture result in flow channelling, which is all the larger as the fracture is more closed and as the correlation length between the facing fracture walls is larger [START_REF] Méheust | Scale effects related to flow in rough fractures[END_REF]. Flow channelling usually coincides with a deviation of the hydraulic conductivity of the fracture from that of the parallel plate of identical mechanical (mean) aperture. Surprisingly, this effect has up to now been largely ignored in most studies dealing with foam flow in rock fractures, with a few striking exceptions, which we present below.

Indeed, relatively few studies on foam flows in fractured media are reported in the literature, and several of them neglect fracture wall roughness, considering Hele-Shaw cell fractures which consist of the space between two parallel plane walls. For example, Yan et al. [START_REF] Yan | Foam sweep in fractures for enhanced oil recovery[END_REF] studied the flow of an aqueous foam in a Hele-Shaw cell fracture at ambient conditions and showed how the gas fraction, mean bubbles size and ratio between the latter and the fracture aperture control foam displacement efficiency. They also showed experimentally, and evidenced theoretically, a diversion mechanism in which the foam tends to flow in thinner fractures rather than thicker fractures within a fractured medium. Osei-Bonsu et al. [START_REF] Osei-Bonsu | Fundamental investigation of foam flow in a liquid-filled Hele-Shaw cell[END_REF] also used a Hele-Shaw cell set-up analogous to a smooth fracture and investigated the link between foam quality, mean bubble size, mean flow behavior and the foam's apparent viscosity; they also measured the recovery factor of a silicone oil. The experimental and theoretical study by Jones et al. [START_REF] Jones | Structure-dependent mobility of a dry aqueous foam flowing along two parallel channels[END_REF] provides detailed insight into how pressure drops develop in such smooth channels when the foam can be considered two-dimensional (a feature observed by Kovscek et al. in natural fractures [START_REF] Kovscek | Engineering, Foam flow through a transparent rough-walled rock fracture[END_REF], see below) due to friction at the Plateau border between the lamella separating bubbles and wall films The effect of surface roughness was also confirmed by Polden et al. [START_REF] Polden | Foam Generation and Flow in Fractures with Different Aperture[END_REF]. Buchgraber et al. [START_REF] Buchgraber | Microvisual investigation of foam flow in ideal fractures: role of fracture aperture and surface roughness[END_REF] went one step further in addressing geometrical complexity, by considering both a constant aperture (Hele-Shaw) fracture, a checkerboard-like fracture with two different apertures of 20 and 40 micrometers, and a fracture possessing a three-level synthetic roughness, to study the effect of aperture variations on foam flow. They found that the geometry of the fracture can greatly impact foam generation mechanisms and foam texture. They visualize and analyzed quantitatively the size distribution of bubbles in situ, but did not measure their velocities.

To our knowledge the first study to address foam flows in a realistic fracture geometry is that by Kovscek et al. [START_REF] Kovscek | Engineering, Foam flow through a transparent rough-walled rock fracture[END_REF], who studied the radial flow of an aqueous foam through a transparent replica of a natural rough-walled fracture at ambient conditions, and measured the pressure drop for various global flow rates and foam qualities [START_REF] Kovscek | Engineering, Foam flow through a transparent rough-walled rock fracture[END_REF]. They could characterize qualitatively the morphology of the bubbles in the fracture plane and in planes perpendicular to the fracture plane, evidencing the two-dimensionality of the foam under these conditions. But they did not obtain a quantitative characterization of the flow from the visualization of the bubbles. Very recently, AlQuaimi and Rosen [START_REF] Alquaimi | Fuels, Foam generation and rheology in a variety of model fractures[END_REF] performed very detailed flow experiments at ambient conditions in synthetic fractures made of roughened glass plates to study foam generation mechanisms and propagation and found that fracture-wall roughness geometry played a major role in controlling the mechanisms of foam generation as well as the foam's mobility. Note however that the geometries which they consider (roughened plates) does not correspond to the known geometry of fractures, i.e., geometries obtained from the fracturing of solid materials (rocks, metals, ceramics); in particular, AlQuaimi and Rossen's geometries exhibit different spatial correlation properties than those of walls resulting from a fracturing process, and are much more akin to two-dimensional porous media than the latter.

Other studies have addressed fractured media at a scale larger than the scale of a single fracture. Haugen et al. [START_REF] Haugen | Miscible and immiscible foam injection for mobility control and EOR in fractured oil-wet carbonate rocks[END_REF] studied foam flow in a carbonate rock cores cut along a plane containing their axis with a saw, at ambient conditions, and found that using pre-generated foam greatly enhanced oil recovery by diverting the injected fluid from fracture to matrix.

Fernø et al. [START_REF] Fernø | Experimental study of foam generation, sweep efficiency, and flow in a fracture network[END_REF] performed foam flow experiments at ambient conditions in fracture networks consisting of fractured marble tiles. They measured the global rheology of the foam and performed a qualitative analysis of bubble geometries some of the fracture. They also found that foam improved sweep efficiency significantly and delayed gas breakthrough. Using a synthetic porous two-dimensional micromodel containing a parallel plate fracture in its middle, Telmadarreie and Trivedi [START_REF] Telmadarreie | New insight on carbonate-heavy-oil recovery: pore-scale mechanisms of post-solvent carbon dioxide foam/polymer-enhanced-foam flooding[END_REF] examined the pore-scale phenomenology of heavy oil displacement by foam in fractured carbonate reservoirs at ambient conditions. Their mostly qualitative analysis revealed that foam not only acts as a mobility enhancement agent but also significantly increase heavy oil recovery from the oil-wet matrix. Gauteplass et al. [START_REF] Gauteplass | Pore-level foam generation and flow for mobility control in fractured systems[END_REF] studied the mechanisms of foam generation in a porous microfluidic micromodel mimicking a sandstone geometry and characterized the role of snap off, evidence that it also occurs at the boundaries between fractures and the porous matrix.

In all the preceding studies addressing foam flow at the fracture scale, the only one to consider a realistic fracture geometry is the seminal paper by Kovscek et al. [START_REF] Kovscek | Engineering, Foam flow through a transparent rough-walled rock fracture[END_REF], but it addresses a radial flow (i.e., a point injection) for which the velocity decreases as the inverse of the distance to the injection point. The prime objective of the present research was to provide insight into the role of aperture spatial variations on the shear rheology of the foam, the foam bubbles size, and foam coarsening in a realistic fracture geometry under parallel flow conditions, with the additional knowledge of the bubble velocity field (following recent studies of foam flows in two-dimensional porous media [START_REF] Géraud | The flow of a foam in a two-dimensional porous medium[END_REF][START_REF] Géraud | Lamella division in a foam flowing through a two-dimensional porous medium: A model fragmentation process[END_REF]), which has been measured in none of the earlier studies on foam flows in fractures. To do so, single-phase foam flow experiments were carried out using a replica of a Vosges sandstone fracture with wellcharacterized aperture map at several flow rates and different foam qualities. The results were compared with those performed by Osei-Bonsu et al. [START_REF] Osei-Bonsu | Fundamental investigation of foam flow in a liquid-filled Hele-Shaw cell[END_REF] in a Hele-Shaw cell of constant aperture equal to the hydraulic aperture of the Vosges fracture.

Materials and Methods

Experimental Setups and Procedure

The experiments were conducted by injecting pre-generated foam in a replica of a Vosges sandstone fracture with a length L of 26.0 cm, a width w of 14.8 cm, a mechanical (mean) aperture h m of 0.86 mm and a hydraulic aperture h of 0.500 ± 0.005 mm. The details about the fabrication of the replica can be consulted elsewhere [START_REF] Isakov | Fluid flow through rough fractures in rocks I: high resolution aperture determinations[END_REF][START_REF] Radilla | Modeling non-Darcian single-and two-phase flow in transparent replicas of rough-walled rock fractures[END_REF]. The aperture map of the fracture replica, which was obtained by Nowamooz et al. [START_REF] Nowamooz | Non-Fickian transport in transparent replicas of rough-walled rock fractures[END_REF] by using an image processing procedure based on the attenuation law of Beer-Lambert, is presented in FIGURE 1 This so-called self-affine behavior is characterized by the exponent of self-affinity (the socalled Hurst exponent ), such that the exponent of the Fourier power spectrum should be 𝐻 [START_REF] Schmittbuhl | Scaling invariance of crack surfaces[END_REF]. Hence here. The properties of apethe rture field of geological

-(2𝐻 + 1) 𝐻 = 0.3
fractures are in large part inherited from those of the facing fracture walls which define them.

In particular, the Fourier spectrum of the aperture field is expected to be self-affine with a Hurst exponent characteristic of the fracture surface topographies (which, for sandstone, is 0.5 [START_REF] Boffa | Roughness of sandstone fracture surfaces: Profilometry and shadow length investigations[END_REF]) at scales smaller than the correlation length between the two fracture walls, and to exhibit a lesser decrease at larger scales due to matching of the two walls at these scales [START_REF] Brown | Correlation between the surfaces of natural rock joints[END_REF][START_REF] Brown | Simple mathematical model of a rough fracture[END_REF]. The measured aperture field is consistent with such behaviour if one considers that the correlation length is on the order of the fracture's length. In addition, an inflection of the spectrum towards a -2 exponent is visible at the largest probed wave numbers (see FIGURE 1(c)), Similarly, the PDF of fracture wall topographies has been measured to be Gaussian [START_REF] Ponson | Failure mechanisms and surface roughness statistics of fractured Fontainebleau sandstone[END_REF], from which it follows that the aperture field should have a PDF that is not very different from a Gaussian. In our data (FIGURE 1 (b)), a log-normal distribution which may be suggested by the asymmetry of the curve has been ruled out by plotting the PDF of log(a), and so the hypothesis of a Gaussian distribution with an additional secondary peak is the most likely. The secondary peak could be due to chemical weathering of the minerals over long times. The results were compared to those obtained by Osei-Bonsu et al. [START_REF] Osei-Bonsu | Fundamental investigation of foam flow in a liquid-filled Hele-Shaw cell[END_REF] using a Hele-Shaw cell with a uniform aperture of 0.5 mm size (equal to the hydraulic aperture of the real fracture) and smooth surface. This allowed the identification of the impact of aperture variations on the foam flow. A customized Hele-Shaw cell with dimensions of 31 × 20 × 0.5 cm was used by the latter authors, and the foam was generated by simultaneously pumping air and aqueous surfactant solution through a foam generator similar to the one used in the present experiments. The pre-generated foam was then delivered to the cell through a tube connected to a metal pipe junction, which was screwed to the inlet of the medium. Pressure transducers were also connected to the inlet and the outlet to measure the pressure drop throughout the Hele-Shaw cell.

At the end of the foam flood experiments in the real fracture, after the injection had been interrupted and once the pressure gradient along the fracture could be considered to have approached zero and foam bubbles to have become stagnant, the evolution over time of the size of the bubbles was recorded in different zones in order to assess the impact of the local aperture on foam coarsening.

Finally, another type of flow cell was used to characterize the bubble size distribution of the foams. Indeed, foams obtained with the same foam generator as in the rough fracture flow experiments were flown through a Hele-Shaw cell, that is, a parallelepipedal tank of thickness much 1 mm much smaller than its two other dimensions, and pictures of the foams were taken.

Image processing & segmentation

A combination of the open source image-processing program ImageJ [START_REF] Rasband | Imagej, us national institutes of health[END_REF], and of either custom-written MATLAB scripts or available MATLAB functions were used to segment and analyze the image sequences captured during the experiments. Two types of analyses were performed.

Firstly, the foams images taken in the 1mm Hele-Shaw cell were treated to measure the bubble size distribution of the foams, for the various foam qualities. The images were first thresholded and transformed into black and white pictures were the bubbles appeared white on a black background (see figure 2(a)), then the equivalent diameter of all connected white regions in the images was measured as from their area . We thus obtained the 𝑑 = (𝐴/𝜋) 1/2 𝐴 probability density functions for the bubbles' apparent diameters. There are power laws (see figure 2(b)), except for the largest foam quality (0.98) for which we obtained an exponential distribution. The foam image appeared very different in the latter case, as compared to the others, so obtaining a different type of bubble size distribution was not surprising. For all images a range of small apparent diameters that were outside of the functional trend and corresponding to connected white regions with shapes very different from the circular shapes, were discarded in the analysis. The bubbles whose apparent diameter is smaller than the HS cell's thickness ( mm) can 𝑒 = 1 be considered to have real diameter equal to the measured apparent diameter, while those with an apparent diameter larger than 1 mm are squeezed between the top and bottom glass plates, so their real diameter (when unconstrained) is different from the measured apparent diameter. We have considered that the bubbles of apparent diameter three times larger than the confining cell thickness had a cylindrical shape, and that their real diameter could therefore be estimated from the apparent diameter as . In the range of

𝑑 real = ( (3/4)𝑒 ( 𝑑 2 )
2 )

1/3 apparent bubble diameters between 1 and 3 mm we have considered a smooth monotonic transition based on a nonlinear combination between these two regimes. The obtained results

do not depend much on the functional form taken for that nonlinear combination. Finally the average diameter for the foam bubbles was computed from the distribution. Given the fact that the distribution is a power law, it was obviously important to perform the measurement on a large population of bubbles.

Secondly, pictures or the foams flowing inside the rough fracture, taken at regular time intervals, were used to obtain maps of velocity vectors of the foam bubbles inside the fracture using PIV (Particle Image Velocimetry) analyses, which estimates local velocity by crosscorrelating sub-regions of the images between two successive time steps. We used package PIVLab [START_REF] Thielicke | PIVlab-towards user-friendly, affordable and accurate digital particle image velocimetry in MATLAB[END_REF] in MATLAB, with a single interrogation size of 128× 18 on the first pass, and 64 × 64 pixels on the second pass. The raw images were first applied a high pass filter of size 15 pixels (see figure 3(a-b)). The horizontal dimension of 2048 pixels corresponded to a length of 260 mm, hence one pixel corresponded to 127 µm in real space, about one order of magnitude smaller than the typical bubble size, but the PIV procedure provided a horizontal resolution of the velocity map of 16 mm. Note that the bubbles whose velocities were measured in this manner belonged to the top layer of bubbles, in contact with the upper wall of the fracture. For some subparts of the image, the filtered image did not contain any significant "dots" to be traced, so no velocity was measured in these regions (see figure 3(bd)). It is, in particular, the case for four circular regions along the center line of the flow cell, where pressure sensors were installed.

Results and discussion

Dependence of the initial mean bubble size on the foam quality

From analyzing images of the foams in the Hele-Shaw cell as explained in section 2.2, we obtain the mean bubble size for the freshly-generated foams, that is, the initial bubble size for the foams that are being injected in the rough fracture. This initial mean bubble size increases monotonically with the foam quality (see figure 4 (a)) but only varies by 25% when the foam quality varies between 80% and 95%. Between 95% and 98% the foam structure changes sharply and consequently the increase in mean bubble size is much more abrupt.

Bubble velocity fields

For Newtonian fluids, it is well known that the presence of in a real fracture of spatiallycorrelated regions of large (respectively, low) apertures results in the creation of tortuous preferential flow paths [START_REF] Brown | Fluid flow through rock joints: the effect of surface roughness[END_REF][START_REF] Méheust | Flow enhancement of a rough fracture[END_REF][START_REF] Méheust | Geometrical heterogeneities and permeability anisotropy of rough fractures[END_REF][START_REF] Brown | Experimental observation of fluid flow channels in a single fracture[END_REF]. This flow channeling impacts the fracture's permeability (or transmissivity), but is also responsible for the anomalous (non-Fickian) transport [START_REF] Nowamooz | Non-Fickian transport in transparent replicas of rough-walled rock fractures[END_REF].

According to the results of Tsang [START_REF] Tsang | The effect of tortuosity on fluid flow through a single fracture[END_REF] for Newtonian fluids, this mechanism can generate pressure drops one or two orders of magnitude higher as compared to the case of a smooth Hele-Shaw cell. In the present experiments, tracing the trajectory of bubbles in the fracture plane shows that the bubbles do not follow an extremely tortuous path that avoids the lowaperture region in the central part of the fracture (see for example figure 3), but their velocity varies significantly along their trajectories, decreasing sharply in this central region. One should remember that such a foam flow is overall incompressible, so if the overall flux for this layer of bubbles varies longitudinal, it must be compensated by an antagonist change in the flux of other bubble layers which we do not measure. In any case, the aperture heterogeneities impact the velocity field significantly, and, therefore, the values of are 𝜇 app expected to be higher in the rough fracture than in a smooth fracture of identical mean aperture.

Dependence of the apparent shear viscosity of the foams on their quality

The relationship between foam quality and for both the real fracture and e Hele-Shaw 𝜇 app cell was calculated from the experimental measurements using Eq. ( 1). The results are presented in FIGURE 4 (b). It can be noticed that is much higher in the real fracture than 𝜇 app in the Hele-Shaw cell for all investigated foam qualities at the considered flow rate, which is consistent with the increase in pressure drop resulting flow channelling which is responsible for changes in hydraulic conductivity. Note that the pressure loss sensitivity to aperture fluctuations is expected to be higher in the case of shear-thinning fluids and yield stress fluids as foams (shear viscosity depends on the local aperture) [START_REF] Roustaei | Non-Darcy effects in fracture flows of a yield stress fluid[END_REF], as compared to Newtonian fluids. and the aperture threshold. The gas flow rate (q g ) is 10 ml/min. It can also be observed from FIGURE 4 (b) that the relationship between the foam's apparent viscosity and its quality in the Vosges fracture is monotonically decreasing whereas this dependency is monotonically increasing in the Hele-Shaw cell. In their work, Osei-Bonsu et al. [START_REF] Osei-Bonsu | Fundamental investigation of foam flow in a liquid-filled Hele-Shaw cell[END_REF] stated that the bubbles of foam with a higher quality offer more resistance against interface deformation and hence exhibit lower mobility in the Hele-Shaw cell. In contrast, is all the larger as the foam quality is smaller in the present case of realistic fracture 𝜇 𝑎𝑝𝑝 geometry, which is in agreement with some previous studies [START_REF] Buchgraber | Microvisual investigation of foam flow in ideal fractures: role of fracture aperture and surface roughness[END_REF][START_REF] Alverez | A unified model for steady-state foam behavior at high and low foam qualities[END_REF][START_REF] Vassenden | Experimental foundation for relative permeability modeling of foam[END_REF]. Note however that foam was mainly generated by snap off in a foam generator. As we increased the flow rate of surfactant to decrease the foam quality, more surfactant solution was available for foam generation, leading to an increased rate of foam generation and a decrease in the size of the mean bubble size of the foam at the inlet of the fracture (i.e., mean initial bubble size), as shown in FIGURE 4 (a). This also explains the decreasing dependence of the apparent viscosity on the foam quality (FIGURE 4 (b)). Note however that in the data of FIGURE μ app 4(c) the foam quality is not varied independently of all other parameters. Part of the dependence of on the foam quality is due to the slight increase in mean velocity μ app observed as the foam quality is decreased, which clearly appears in FIGURE 4 (c).

Visual observations (which will be discussed in more details in section 3.4) revealed that foam bubbles were stagnant in some parts of the fracture, which means (not surprisingly) that the foam behaved as a yield stress fluid. These observations also showed, in agreement with previous studies [START_REF] Princen | Rheology of foams and highly concentrated emulsions: I. Elastic properties and yield stress of a cylindrical model system[END_REF][START_REF] Höhler | Rheology of liquid foam[END_REF], that the extension of the stagnant area in the rough fracture becomes smaller as the foam quality was decreased and smaller foam bubbles were generated. The apertures encountered in these areas where bubbles were stagnant were smaller than a threshold value. This threshold value was observed to vary with the foam quality, as shown in FIGURE 4(d), due to the dependence of the size of the bubble on the foam quality (FIGURE 4(a)). A remarkable feature of the flow patterns in the real fracture is that bubbles were stagnant in the lowest aperture regions. As expected, the foam behaved as a yield stress fluid which does not flow when the applied shear stress is below a threshold value. In the present case, the local shear stress in the low-aperture area fell below the yield stress as the pressure gradient was not high enough at the lowest flow rates. The required pressure gradient to exceed this threshold is obtained by increasing the flow rate. This phenomenon can be seen by comparing the velocity map at 10 and 20 ml/min in FIGURE 5. As the flow rate increased, a larger area of the fracture plane was occupied by flowing foam bubbles.

Shear rate distribution in the rough fracture

Shear rate distribution in the fracture plane

Given that the foam is a shear thinning fluid, the variation in velocity leads to different viscosity values in different regions of the fracture. The equivalent Newtonian wall shear rate was calculated from the aperture and velocity magnitude maps using the following 𝛾 𝑁 expression [START_REF] Pipe | High shear rate viscometry[END_REF].

(2)

𝛾 𝑁 = 6𝑣 𝑎
In the preceding equation, v is the local velocity magnitude, i.e. the average value of the velocity profile between both surfaces of the fracture. Here we assume that the local velocity magnitude which we measure and is represented in FIGURE 5 can be assimilated to v in Eq.

(2). In the same equation, a is the local aperture value as represented in the aperture map of Figure 1. Eq. ( 2) applies to the 2-D flow of a constant viscosity liquid between two infinite parallel planes (the plane Poiseuille configuration). However, for flows of liquids with a shear-rate-dependent viscosity, the calculation of the shear rate is more complex as the velocity profile is no longer parabolic [START_REF] Macosko | Rheology: principles, measurements, and applications[END_REF]. In this case, the true wall shear rate can be found using the Weissenberg-Rabinowitsch-Mooney equation [START_REF] Macosko | Rheology: principles, measurements, and applications[END_REF]:

, 𝛾 = 𝛾 𝑁 3 [ 2 + 𝑑(𝑙𝑛𝛾 𝑁 ) 𝑑(𝑙𝑛𝜏) ] (3)
with the wall shear stress. It must be noted that the rough approximation has been 𝜏 𝛾 ≈ 𝛾 𝑁 made in the present analysis. Improved accuracy can be achieved by calculating , which

𝑑(𝑙𝑛𝛾 𝑁 ) 𝑑(𝑙𝑛𝜏)
is needed as input in Eq. ( 3). However, the experimental pressure distribution throughout the fracture would be required for the calculation of , which was not available in the present 𝜏 experiments.

The results in terms of wall shear rate distribution are shown in FIGURE 5, for foams of quality 85% and flow rates of 10 and 20 ml/min. This figure shows that the shear rate in the rough fracture is strongly impacted by the spatial variations of the aperture field, as expected from its definition. In contrast, it would be uniform within the Hele-Shaw cell of [START_REF] Osei-Bonsu | Fundamental investigation of foam flow in a liquid-filled Hele-Shaw cell[END_REF]. The zero shear rate regions are simply the regions where the velocity of the foam is zero, which exist due to the fact that the foam is a yield stress fluid. The extension of these regions of stagnant bubbles became smaller as the flow rate was increased (see Figure 5).

The steady-state shear flow of bulk foams has been proved to be well described by the Herschel-Bulkley law [START_REF] Herschel | Konsistenzmessungen von gummi-benzollösungen[END_REF] in previous works [START_REF] Reidenbach | Advertising and male nudity: An experimental investigation[END_REF][START_REF] Harris | Rheological properties of low-gel-loading borate fracture gels[END_REF][START_REF] Bonilla | Experimental investigation on the rheology of foams[END_REF]. This empirical law can be written as follows:

{ 𝜏 = 𝜏 0 + 𝑎𝛾 𝑛 for 𝜏 > 𝜏 0 𝛾 = 0 for 𝜏 ≤ 𝜏 0 (4)
Where is the yield stress, is the consistency and n is the flow index of the fluid. In the 𝜏 0 𝑎 case of shear-thinning yield stress fluids, n is inferior to unity. For bulk foams it is usually in the range (see for example [START_REF] Denkov | Wall slip and viscous dissipation in sheared foams: Effect of surface mobility[END_REF]). The three parameters are generally calculated by [0.2;0.4] fitting the data obtained by measuring the shear rate as a function of the applied shear stress 𝛾 τ using a rheometer, which is particularly challenging in the case of foams. For a foam flowing through a porous medium, the rheology is also expected to follow a similar Herschel-Bulkley rheology, but with values of the yield stress, consistency and yield stress that are different from those exhibited in bulk foam flow. Our present study provides hints towards a method to measure the in situ rheology of the foam flowing through the porous medium.

The apparent shear rate of the complex fluid flowing through the real fracture can be 𝛾 𝑎𝑝𝑝 obtained from the ratio of the characteristic pore velocity q (obtained by dividing the total flow rate by the cross-sectional area) and the characteristic microscopic length usually taken as : 𝐾 𝛼 the bulk rheology of the fluid and the porous media [START_REF] Chauveteau | Rodlike polymer solution flow through fine pores: influence of pore size on rheological behavior[END_REF][START_REF] Sorbie | The rheology of pseudoplastic fluids in porous media using network modeling[END_REF][START_REF] Lopez | Predictive network modeling of single-phase non-Newtonian flow in porous media[END_REF][START_REF] Comba | Rheological characterization of xanthan suspensions of nanoscale iron for injection in porous media[END_REF].

is commonly used as an 𝛾 𝑎𝑝𝑝 input in the rheological model (Eq. 4) in order to predict the viscosity of the fluid in the porous medium. Nevertheless, , and n are unknown in the case of foams, so the preceding 𝜏 0 𝑎 approach is hardly applicable. Also, it has been traditionally assumed that is independent of 𝛼 q. However, Rodríguez de Castro and Radilla [START_REF] Rodríguez De Castro | Flow of yield stress and Carreau fluids through rough-walled rock fractures: Prediction and experiments[END_REF] recently showed, using flow experiments of polymer solutions in rough fractures, that this assumption is no longer valid in the presence of a yield stress.

The average wall shear rate throughout the media, obtained from averaging the shear rate 𝛾 𝑎𝑣 maps such as those shown in FIGURE 5 is plotted as a function of the foam quality for a flow rate of 10 ml/min in Figure 6(a). FIGURE 6(a) shows that the average shear rate is not monotonic as a function of the foam quality in the real fracture. It is reminded here that foam quality was varied by setting a different liquid flow rate while keeping the gas flow rate constant. Therefore, decreasing foam quality leads to increased rate of foam generation with finer texture (smaller bubbles). According to Princen [START_REF] Princen | Rheology of foams and highly concentrated emulsions: I. Elastic properties and yield stress of a cylindrical model system[END_REF] and Hohler and Cohen-Addad [START_REF] Höhler | Rheology of liquid foam[END_REF], the yield stress is inversely proportional to the size of the bubbles. Finer foam bubbles at lower foam quality have lower yield stress that lead to the flow of foam bubbles to a wider area of the fracture (FIGURE 4 (b-d)), which tends to increase the average shear rate. On the contrary, decreasing the foam quality means decreasing the average flow velocity (see Figure 4(c)) to some extent, which tends to promote a lower mean shear rate. These two competing effects may explain the non-monotonicity of the plot in Figure 6(a).

The average shear rate was also computed for the series of experiments at constant foam quality 85% and various flow rates ranging between 10 and 60 ml/min, but the acquisition frequency of 1 Hz allowed obtaining meaningful velocity fields only for the data recorded at 10 and 20 ml/min (at 11.5 and 22.6 m/s, respectively), with an increase when increasing the flow rate, as expected. In contrast, the dependence of the apparent viscosity could be obtained from Eq. ( 1)) for all flow rates. Figure 6(b) shows apparent viscosity decreases with flow rate in both fracture and Hele-Shaw cell as expected due to shear thinning behavior of foam. It is shown in Figure 6(c) and seems consistent with a power law behavior, of exponent -0.41, to be compared to the exponent -0.27 for the Hele-Shaw cell data of Osei-Bonsu et al. Since only two measurements of the average shear rate were available, we could not check that the dependence of the apparent viscosity on the average shear rate was consistent with a power law behavior, but we propose that a method based on such "apparent rheograms" could be used to measure the in situ rheology of the foam in the porous medium. Similarly, the shift factor ( assuming in Eq. ( 5)) could be computed for 𝛼 = 𝛾 𝑎𝑣 𝐾 /𝑞 𝛾 𝑎𝑣 ~𝛾𝑎𝑝𝑝 the two flow rates 10 and 20 ml/min, to 1.30 and 1.26, respectively.

Correlations between foam coarsening and aperture distribution

The nature of foams is dynamic; their structure changes with time due to foam coalescence by liquid drainage, film rupture and gas diffusion across the films separating neighboring bubbles [START_REF] Höhler | Rheology of liquid foam[END_REF]. It is interesting to observe how the local fracture aperture influences the coarsening dynamics of the foam after the flow has been stopped. To this end, closed views of the foam's structure were recorded in sub-windows of the entire fracture plane positioned at different positions in that plane. All these sub-windows corresponded to different mean apertures (computed over the sub-window). The evolution of the mean bubble sizes was monitored for a total duration of 15 hours. FIGURE 7(a) shows pictures taken for three subwindows for which the mean aperture was respectively 0.2, 0.5 and 1.5 mm, at three different times after the flow has been stopped: 0 s, 3000 s, and 15 h. The coarsening in time of the foam is clearly visible in all image, and appears to be all the faster as the mean aperture is larger. FIGURE 7(b) shows the corresponding dependence of the mean bubble size on 〈𝑑〉 sw the mean aperture , at times 0 s and 15 h after the flow has been stopped. The gas 〈𝑎〉 sw diffusion rate through the films is expected to be higher at earlier times as smaller foam bubbles with higher internal pressure exist in the system. Gas diffuses from smaller to larger bubbles with a rate that is proportional to the pressure difference between them, eventually leading to the vanishing of the smallest bubbles and continuous coarsening of the foam. The mean bubble size is indeed observed increase in time, for a given , as expected, 〈𝑑〉 sw 〈𝑎〉 sw but as we only have measurements at two different times we cannot quantify the changes in coarsening rates. More interestingly, at a given time the mean bubble size appears to be highly correlated to the mean aperture in the sub-windows, following a monotonically increasing dependence.

It is interesting to note that the mean aperture for which the plots are above the 〈𝑑〉 sw = 〈𝑎〉 sw line correspond to regions in which the mean horizontal size of the bubbles is larger than the aperture that confines them vertically. Consequently, it is very likely that there is only one layer of bubbles in these regions of the fracture planes. In contrast, as the mean aperture is larger than the mean horizontal bubble size, the bubbles are expected to be pseudo-spherical as they are not confined between the two fracture walls. Therefore, more than one layer of bubbles can exist in these regions. The plot obtained for the initial time (after the flow has been stopped) is particularly interesting as it shows the structural state of the foam resulting from its flow in the fracture. This means that smaller bubbles flow in smaller aperture regions, while larger bubbles flow in larger aperture regions. Furthermore, only bubbles present in low aperture regions (or mean aperture smaller than 0.5 mm) are confined vertically (despite being smaller than in regions of larger aperture). On the contrary, at t = 15 h after the flow has been stopped, the local mean bubble size is larger than the local mean aperture in all sub-windows, and so the foam can be expected to be confined vertically in all of them. In other words, at t = 15h the foam is two-dimensional in the entire fracture plane (i.e., there is a monolayer of bubbles everywhere). As Cohen-Addad [START_REF] Cohen-Addad | Viscoelastic response of a coarsening foam[END_REF] showed, the foam's apparent viscosity decreases during foam coarsening process. Accordingly, we can expect a deterioration of the sweep efficiency when using coarsened foam as displacing fluid, given the decrease in .

𝜇 app

Conclusions

As for the creeping flow of a Newtonian fluid [START_REF] Méheust | Geometrical heterogeneities and permeability anisotropy of rough fractures[END_REF][START_REF] Brown | Fluid flow through rock joints: the effect of surface roughness[END_REF], spatial variations of the aperture field of a rough fracture with a realistic geometry strongly impacts the flow of an aqueous foam.

To the best of our knowledge, all previous experimental studies addressing foam flow in fractured media had only been conducted on fractured media manufactured by glass plates [START_REF] Yan | Foam sweep in fractures for enhanced oil recovery[END_REF][START_REF] Osei-Bonsu | Fundamental investigation of foam flow in a liquid-filled Hele-Shaw cell[END_REF][START_REF] Buchgraber | Microvisual investigation of foam flow in ideal fractures: role of fracture aperture and surface roughness[END_REF][START_REF] Perez | Foam Injection into horizontal cells with channels of different openings[END_REF], with unrealistic fracture geometries, or real fracture media but with no quantitative characterization of the aperture field [START_REF] Haugen | Miscible and immiscible foam injection for mobility control and EOR in fractured oil-wet carbonate rocks[END_REF][START_REF] Kovscek | Foam flow through a transparent rough-walled rock fracture[END_REF][START_REF] Brattekas | Fracture mobility control by polymer gel-integrated EOR in fractured, oil-wet carbonate rocks[END_REF][START_REF] Ferno | Experimental study of foam generation, sweep efficiency and flow in a fracture network[END_REF]. In additional, though several of these previous studies featured local visualization of bubble morphologies [26][40, 41, 44] [ 45, 46][48], and in some cases the measurement of bubble size distributions [START_REF] Yan | Foam sweep in fractures for enhanced oil recovery[END_REF][START_REF] Buchgraber | Microvisual investigation of foam flow in ideal fractures: role of fracture aperture and surface roughness[END_REF], none of them reported bubble velocity measurements nor tried to link local bubble velocities to the in situ rheology of the foam. This study set out to provide an understanding of the effects of spatial aperture variations on the flow patterns, shear rheology, and bubble morphology during the flow of foams through a realistic fracture geometry, that is, a geometry that possesses the right spatial correlation properties which are known to be typical of fracture surfaces and the gap between them [START_REF] Brown | Simple mathematical model of a rough fracture[END_REF][START_REF] Brown | Correlation between the surfaces of natural rock joints[END_REF]. To do so, a comprehensive series of single-phase flow experiments were conducted in a replica of natural rough-walled fracture (Vosges sandstone), and confronted to similar measurement performed in a Hele-Shaw cell fracture of identical hydraulic aperture. Based on these measurements, the following observations and conclusions can be drawn:

 The aperture variation strongly affects the apparent viscosity of foam. In particular, the existence of spatially-correlated low-and high-aperture areas increases flow tortuosity and the velocity contrast between preferential flow channels and stagnant zones, as is the case for the flow of a Newtonian fluid. By modifying the velocity field and, consequently, the shear rate distribution, it controls the in situ rheology of the foam.

 Foam is a yield stress fluid and exhibits shear-thinning behavior when flowing through a porous medium, but the rheology is not necessarily the same as that of the bulk foam. Here we could not obtain a complete "in situ" rheogram for the 85 % foam from foam bubble velocity measurements, but we propose to use this method to measure the in situ rheology of the foam. In this case, the measurements that could be made point to a power law of exponent -0.41.

 The foam's bubble size inside the fracture is correlated to the local fracture aperture, and is all the larger as the local aperture is larger. The transition from bulk foam to a vertically-confined two-dimensional foam occurs in the fracture plane along the equiaperture lines a=0.5 mm in our setup.

 The dynamics of foam coarsening inside the fracture is correlated to the local aperture and is all the faster as the local aperture is larger.

A replica of a Vosges sandstone fracture with a given Hurst exponent was considered in the present study to investigate the effect of spatial distribution on apparent shear viscosity and bubble morphology. In future studies, it would be interesting to consider similar geometries with various Hurst (H) exponent values as well as, possibly, different matching (i.e., correlation) scales between the two fractures walls, and correlate the flow properties of the foam to theses stochastic geometrical parameters.

  (a). The minimum aperture is 0 (closed regions) and the maximum aperture 2.294 mm. This figure shows that the smallest apertures are located in the center of the model, while the largest apertures are close to the inlet and the outlet. The probability density function (PDF) of the fracture apertures, shown in FIGURE1(b), is peaked around 0.71 mm but not symmetrical, as a second smaller peak is seen around 1.13 mm. The power spectrum of one-dimensional profiles of the aperture map parallel to (FIGURE 1(c)) exhibits a power law behavior with a x global exponent -1.6, which indicates long-range spatial correlations along the aperture field.

FIGURE 1 .

 1 FIGURE 1. (a) Aperture map of the fracture used in the present experiments (data obtained

FIGURE 2 .

 2 FIGURE 2. (a) Image of foam of quality 85% in the Hele-Shaw cell (with 12535 bubbles).

FIGURE 3 .

 3 FIGURE 3. (a) Raw image of the foam flowing through the fracture (flow rate 10 ml/min,

FIGURE 4 .

 4 FIGURE 4. (a) Dependence of the mean initial size of the generated bubble on the foam

FIGURE 5 .

 5 FIGURE 5. Maps of velocity magnitude and wall shear rate at 85% foam quality for the

  and is an empirical shift factor known to be a function of both 𝐾 = ℎ 2 12

FIGURE 6 .

 6 FIGURE 6. (a) Dependence of the average shear rate on the foam quality at a 10 ml/min gas

FIGURE 7 .

 7 FIGURE 7. (a) Bubbles observed in different sub-windows of the aperture plane, at

Acknowledgment

We would like to acknowledge the UK Engineering and Physical Sciences Research Council (EPSRC) for providing the Ph.D. studentship for Mohammad Javad Shojaei. We thank Pr.

Giovanni Radilla for providing the equipment used in the experiments. We also thank Frédéric Bastien for his assistance with the experimental setup.


est destinée au dépôt

Thermally Stable Ohmic Contact to p-Type 4H-SiC

Based on Ti3SiC2 Phase