Lionel Mercury

Emmanuel De Bilbao

Patrick Simon

Hugues Raimbourg

Isabelle Bergonzi

Claudie Hulin

Aurélien Canizares

Kirill I Shmulovich

Patrick Simon -Cemhti

Quartz Stressing and Fracturing by Pore Pressure Dropping Down to Negative Pressure

Keywords: tensile strength, pore rock damage, storage safety, Raman scattering, quartz weakening, fluid inclusions

In water-bearing porous rocks, pore pressure variations play a major role in deformation, through dissolution-precipitation and fracturing processes. An often-overlooked variation where pressure falls to negative pressure or tension can operate whenever aquifer formations dry out, for instance, in deep storage (nuclear or industrial wastes, long-term CO 2 mitigation, short-term energetic resources, etc.). This can generate capillary tension within the aquifers. This study investigates the mechanical effect of such in-pore tension in the surrounding crystal field, through laboratory experiments at the one-pore scale. Microthermometric procedures were carried out on synthetic fluid inclusions to generate large tensile stress and were combined with Raman microspectrometry to visualize the resulting stress fields in the host quartz. For comparison, we numerically modeled the stress field by linear elasticity theory. The experiments demonstrate that significant damage is produced in crystalline materials by the pore tension. Despite the induced stress measured by micro-Raman spectrometry to remain moderate, it is able to fracture the quartz. The volume of the cavity is a prominent controlling parameter for the stress amplitude. The crystalline heterogeneities of the solid are another major parameter for localizing the mean weak stress and accumulating overstress. Our results call for bringing pore-scale micromechanics into the safety assessment of the geological storage of various wastes inside depleted aquifers. They also show the magnifying effect of heterogeneities on propagating stress and localizing it along certain directions, promoting the final failure of water-bearing minerals, rocks, or pore networks.

INTRODUCTION

Modern societies have growing interest in water-depleted aquifers for various storage operations, which are often accompanied by enhanced drying and increasing water tension (liquid negative pressure) in the liquid infilling the pores. Water tension in pores does induce a stress in the solid lattice around the fluid-bearing cavity, which is largely overlooked while possibly conducive to crystal failure. Even in the subcritical cracking regime, low stresses, possibly amplified at the tip of microcracks because of pre-existing flaws in the rock, and slow processes, persisting in the long term, can result in rock failure (e.g., refs. 1,2 ).

The micromechanics of such situations can be treated similarly to the drying-driven damage occurring in rocks at variable humidities (e.g., refs. 3,4 ). The theory involves the role of variable water saturation, of capillary bridging forces crossing the micro/nanopores, and of disjoining forces that hinder adsorption of water in nanopores, all processes that cause strain softening and microcracking. In all these studies, it is difficult to distinguish at the relevant scale the mechanical role of field stress created by pore tension from stress corrosion sequences of various fluid injections (air, liquid, and oil).

The water invasion in dried rocks is known to induce weakening effects without a clear viewpoint on the active mechanisms (e.g., refs. [START_REF] Turuntaev | Formation fracturing by pore pressure drop (laboratory study)[END_REF][START_REF] David | Mechanical instability induced by water weakening in laboratory fluid injection tests[END_REF] ). Fourth, salt damage to cultural heritage stones is mostly interpreted as the role of crystallization pressure but is accompanied by significant desaturation and capillary effects (e.g., refs. [START_REF] Li | Capillary tension theory for prediction of early 1001 autogenous shrinkage of self-consolidating concrete[END_REF][START_REF] Ro̷ Yne | Crack propagation driven by crystal growth[END_REF] ).

The key objective of this contribution is therefore to address drying-driven microcracking mechanisms, which may downgrade the material's quality and impact the physical properties (transmissivity and safety storage) of the damaged rock. Two experiments were designed with fluid inclusions synthesized in monocrystalline quartz and filled with aqueous solutions. As such, they stand for pores in crystalline and cemented materials but cannot be directly used to draw conclusions on natural fluid inclusions that contain less structural water and lower dislocation densities than synthetic fluid inclusions (SFIs). The occluded liquid was put under variable high tensions with various time frames, according to a microthermometric procedure (closed univariant system) (e.g., refs [START_REF] Green | Water and solutions at negative pressure: Raman spectroscopic study to -80 Megapascals[END_REF][START_REF] Zheng | Liquids at large negative pressures: water at the homogeneous nucleation limit[END_REF][START_REF] Alvarenga | Elastic properties of water under negative pressures[END_REF][START_REF] Shmulovich | Experimental superheating of water and aqueous solutions[END_REF][START_REF] Fall | The effect of fluid inclusion size on determination of homogenization temperature and density of liquid-rich aqueous inclusions[END_REF] ). The For each inclusion, two stereoplots show the orientation of the [0001] and [112̅ 0] axes within the host quartz, located perpendicular and parallel to the polishing plane, respectively. The redundant coordinate system (a 1 , a 2 , a 3 , c) is used (hexagonal lattice).

ACS Earth and Space Chemistry

Raman shift of the quartz was mapped and shown to vary as a function of the in-pore tension, therefore visualizing the internal stress fields in the host crystal.

Two levels of damage occurred during the course of these experiments: optically invisible damage at the microscale and visible fracturing of the host quartz. The experimental investigations were compared to a numerical simulation based on isotropic linear elasticity. The distance of influence of the in-pore tension inside the host quartz appears strikingly larger than expectable, while the simulated and measured distributions of stress in the crystal match.

MATERIALS AND METHODS

2.1. Samples. SFIs have been synthetized in two monocrystalline pieces of hydrothermal quartz. After thermal fracturing, the quartz pieces were healed at high pressure and temperature (P,T) by SiO 2 overgrowth, in the presence of an aqueous solution of known composition. [START_REF] Shmulovich | Experimental superheating of water and aqueous solutions[END_REF] One quartz piece was healed in the presence of a CsCl 13.6 M (69.6 wt %) brine, at 630 MPa and 480 °C for 20 days. Under ambient conditions, CsCl inclusions contain either two phases (a bubble of gas and a liquid phase) or a monophasic liquid-only f1 phase (Figure 1). This trapped solution is known to sustain the most extreme tensile strength ever recorded in aqueous fluids [START_REF] Shmulovich | Experimental superheating of water and aqueous solutions[END_REF] and therefore will serve to install super-high tensile stress. In addition, the decreasing water activity associated with the CsCl electrolyte (0.57 from Pitzer calculations at 25°C [START_REF] Holmes | Thermodynamic properties of aqueous solutions of the alkali metal chlorides to 250[END_REF] ) may contribute to water weakening by enhancing solubility, especially at the stress concentration zones.

The other quartz piece was hydrothermally treated at 500 °C, 600 MPa for 45 days, in the presence of a 0.2 M NaOH dilute solution. Under ambient conditions, these NaOH inclusions also contain either two phasesa bubble of gas and a liquid phaseor a monophasic liquid-only phase (Figure 1). The dilute NaOH solution trapped under high-density conditions is prone to install liquid tension at room temperature (RT) [START_REF] Shmulovich | Experimental superheating of water and aqueous solutions[END_REF] and so serves to exhibit long-term tensile stress.

These two electrolytes make the quartz more soluble, especially under the synthesis conditions (0.63 GPa, 480 °C, and 20 days for the first; 0.6 GPa, 500 °C, and 45 days, for the second), so they may produce inclusions with very good inner surfaces well prone to hosting superheating processes. According to previous datasets, [START_REF] Shmulovich | Experimental superheating of water and aqueous solutions[END_REF] and as mentioned above, CsCl solutions can sustain the most extreme tensile strength ever recorded in aqueous fluids, while dilute NaOH solutions are efficiently superheated at RT or below (high-density samples). Additionally, the role of both electrolytes during these experiments should have enforced the tensile stress by a chemical "weakening" effect. However, considering the slow kinetics of the dissolution-precipitation process under the thermometric cycle conditions, this effect is certainly of secondary importance.

Electron back-scattered diffraction (EBSD) was used to determine the crystallographic orientation of the quartz samples around the fluid inclusions, after the microthermometric studies, to prevent any electron-driven defects. To do so, thin sections were chemically polished with a colloidal silica suspension (0.04 μm colloidal silica suspension, Struers) and then carbon-coated to prevent charging effects. All sections were tilted by 70°to the electron beam to produce clear diffraction patterns. Data were collected using an EDAX PEGASUS EDS/EBSD system and processed with OIM DC [START_REF] Wagner | The IAPWS Formulation 1995 for the Thermodynamic Properties of Ordinary Water Substance for General and Scientific Use[END_REF] Each isochore is labeled with the molar volume in cm 3 •mol -1 . The CsCl inclusion micrographs illustrate the fluid content present in the cavity at each characteristic step.

ACS Earth and Space Chemistry

Spectra were treated with Renishaw Wire 3 software. The 261 Raman spectra were not baseline-corrected but directly fitted 262 with a Gauss-Lorentz (G/L hereafter) profile, to involve even 263 the tiniest distortions in the analysis. Very small differences between mono-and biphasic conditions could thus be probed because the spectra were acquired with strictly identical conditions (optical configuration and temperature), apart from the pressure inside the inclusion. No Bose-Einstein thermal correction was performed 43 because Raman shifts were compared at identical temperatures under two different pressure conditions (monophasic and biphasic regimes). 39 2.4. Method 3: Super-Stressing the Container toward Failure. Two methods were used to establish whether the stress induced in quartz by the in-pore tension was able to provoke the quartz failure. First, the characteristics of the CsCl fluid inclusion allowed us to subject it to high tension cyclically. The procedure was simply to put the inclusion at various temperatures of measurements (T MEAS ) corresponding to different high tensions and leave the inclusion in this monophasic metastable state until cavitation occurs, bringing back the inclusion to the LV equilibrium state. Then, the inclusion was homogenized and again cooled down to a certain T MEAS pursuing the cycle. From a mechanical point of view, these repeated stress conditions can damage the solid depending on the lattice's defaults, the intensity of the stress, and the number of cycles.

Second, the NaOH fluid inclusions allowed us to install high-tension conditions at RT, so the liquid tension could be readily maintained over a long time. Therefore, stress is accumulated within the sample over time, as long as bubbles are absent (easily verified by optical observations) in the superheated inclusions: this was a permanent or "enduring" tension experiment, corresponding to purely static stress conditions.

Method 4: Simulating the Stress in the Host

Crystal as a Function of in-Pore Fluid Tension. Numerical simulations of the experiments were performed by finite element analysis (FEA hereafter) using Comsol Multiphysics software for a deeper understanding of the stress change observed by Raman spectroscopy. In the first approach, the sample was modeled as a 2D axisymmetric solid containing an oblate spheroidal inclusion. [START_REF] Gregora | Raman Study of AlPO4 Berlinite at the Alpha-Beta 1110 Transition[END_REF] The stress field around a cavity has been studied many times. In particular, Eshelby 45 studied how an applied stress, uniform at large distances, is disturbed by a cavity. From this early work, many analytical or closedform analytical solutions have been proposed for spherical or spheroidal inclusions and considering isotropic or anisotropic materials. [START_REF] Eshelby | The determination of the elastic 1115 field of an ellipsoidal inclusion, and related problems[END_REF][47][48][49][50][51][52][START_REF] Meng | Extending Esh3D Code to Solve Interacting Eshelby's 1133 Inhomogeneity Problems[END_REF][START_REF] Meng | Evaluation of the Eshelby 1135 Solution for the Ellipsoidal Inclusion and Heterogeneity[END_REF] The problem becomes more complicated for a spheroidal cavity hosted in an anisotropic solid. On the other hand, FEA has proved its reliability for solving mechanical problems especially in the case of linear elastic problems but also in the case of geometric nonlinearity and with nonlinear isotropic elastic materials. 55 Actually, FEA can be applied whatever the geometry of the inclusion.

In a first approach, the sample was modeled as a 2D axisymmetric solid containing an oblate spheroidal inclusion, [START_REF] Gregora | Raman Study of AlPO4 Berlinite at the Alpha-Beta 1110 Transition[END_REF] and the finite-element-based results were compared with theoretical results obtained using the code recently developed by Meng et al. [51][52][START_REF] Meng | Extending Esh3D Code to Solve Interacting Eshelby's 1133 Inhomogeneity Problems[END_REF] which gives quasianalytical strain/stress fields following Eshelby's solution, with the elliptic integral approximated through a numerical routine. The purpose of this comparison was to determine the size of the solid surrounding the cavity to avoid boundary effects and define the meshing strategy to avoid numerical artifacts or errors and to better understand the stress field building around the inclusion. An isotropic elastic material was defined in the solid part with Young's modulus of 100 GPa and Poisson's coefficient of ACS Earth and Space Chemistry 0.06. 56 These elastic properties correspond to equivalent quasiisotropic quartz obtained by Voigt-Reuss-Hill averaging. 56 The inclusion was defined with null elastic constants, and a negative pressure P FI (e.g., P FI = -1 MPa) was applied on the cavity wall. For the theoretical solution, the superposition principle was applied as the code because it does not directly f3 apply pressure on the cavity wall (Figure 3, top, left). Unstructured triangular quadratic elements were applied for the FEA, and zero displacement was applied to the symmetry plane (horizontal) perpendicular to the revolution axis (vertical). The stresses σ ii /P FI lying along the radial direction were compared (Figure 3). For the sake of simplicity, the stress sign convention commonly used in mechanics has been changed in this paper, and a positive stress refers to a compressive state. The pressure has also been changed P T = Tr(σ ̿ )/3 instead of P T = -Tr(σ ̿ )/3 to keep its usual meaning, namely, a positive value refers to a compressive state. The stresses were normalized by dividing each of them by P FI ( σ σ *= P / ii ii FI ), and the radial coordinate was divided by the halflength of the ellipse a (r* = r/a).

We considered three shapes defined by the semiaxis ratio b/ a (hereafter called shape factor), where a and b are the semiaxes in the radial and axial directions, respectively. The third semiaxis c is equal to the radial axis in the oblate spheroid. Simulated results are in very good agreement with theoretical results for all the stresses, proving that the solid model size was large enough and the meshing strategy was suitable. The FEA model was validated. The results show clearly the stress concentration effect induced by the cavity and the dependence of the concentration amplitude on the shape applied on the cavity wall (negative pressure PFI) (Figure 4). Stresses that were noted by FEA were calculated using Comsol Multiphysics software. For all the calculations, Young's modulus and Poisson's coefficient were 100 GPa and 0.06, respectively, according to ref. 56 ACS Earth and Space Chemistry 386 anisotropic elastic moduli as for the previous simulation were 387 used. First, the Raman maps all around a CsCl inclusion were 396 recorded at 120 °C while the internal pressure in the liquid was 397 either 0.12 MPa, corresponding to biphasic liquid-vapor 398 assemblage, or -194 MPa, corresponding to the metastable 399 monophasic state with liquid under tension. Each mapping 400 took 20 h. The chosen temperature of 120 °C was sufficiently 401 above the T IN to avoid any nucleation during the Raman 402 measurement time. Indeed, the metastable lifetime depends 403 exponentially upon the distance to T IN.

EXPERIMENTAL RESULTS AND INTERPRETATION

58-60

Two A 1 Raman modes of SiO 2 with highest wavenumbers were studied: the mode at 464 cm -1 and the mode at 206 cm -1 . The signal was mapped in the quartz around the inclusion (width 70 μm × height 110 μm) that contained the fluid. A G/L profile was fitted with the recorded spectrum at each point, and the characteristic wavenumber of the G/L profile was used to build the map. The map obtained with the 464 cm -1 mode in the stable state (biphasic liquid-vapor assemblage) shows actual wavenumbers ranging from 462.5 to f7 463.1 cm -1 with a mean value of 462.8 cm -1 (Figure 7). Black To extract the mean Raman shift provoked by the in-pore 468 tension in the actual sample over the experimental area, the 469 whole data set collected with each map was redrawn as the 470 histogram of the number of occurrences in each wavenumber 471 range, for the two Raman mode maps (464 and 206 cm -1 ) f8 472 under the two stressed conditions (Figure 8). This shows how 473 the wavenumbers measured vary within the experimental area 474 for the two Raman modes and also shows how regular the 475 datasets are.

ACS Earth and Space Chemistry

476

The liquid-vapor series gives the Raman modes at 477 saturation pressure: the Gaussian maxima are at 462.8 and 478 199.9 cm -1 . These measurements were double-checked to test 479 the reproducibility of the recorded values (blue and green 480 histograms), and both datasets match well. The match is 481 strictly perfect for the high wavenumber mode (same main 482 wavenumber and same width) while the low wavenumber A 1 483 mode gave a 0.015 cm -1 change for wavenumber (showing the 484 high stability of the experiment) with the total width varying 485 by a maximum of 0.1 cm -1 . The main information we take 486 from these graphs is the Raman gap that exists between the 487 two 120 °C datasets at different pressures. Each mode (464 488 and 206 cm -1 ) has a specific gap, quite moderate with respect 489 to the precision of the spectrometer, but largely beyond any 490 instrumental effect. Note that here the key factor of the 491 spectrometer is its stability, not its spectral resolution. Specific 492 care was taken to improve stability by avoiding thermal 493 fluctuations and then any human presence in the experimental 494 room. Now, using the calibration of the Raman shift with 495 pressure, 42 it is possible to estimate the level of strain involved. 496 The two modes agreed remarkably well for the resulting stress: 497 20 MPa ± 2.5 MPa (Figure 8). The small difference between 498 the values (22.2 MPa for 464 cm -1 mode and 17.6 for 206 499 cm -1 mode) can be related to the wavenumber measurement 500 uncertainty, and by the application of only one parameter on a 501 noncubic structure, which is, moreover, not subjected to a 502 hydrostatic pressure but to a complex stress field.

503

The Raman shift difference of the 464 cm -1 mode between 504 the biphasic and monophasic modes was plotted as a map f9 505 (Figure 9). This plot is a spatial representation of the Raman wavenumber change, that is to say the stress distribution induced by the in-pore pressure change. The difference is positive or zero, except for very rare occurrences. In addition, the average wavenumber shift and standard deviation were evaluated in four zones surrounding the inclusion. The area above the inclusion in the map is called "North," and so on.

The statistics are very close in the three zones West, North, and East (see Figure 9) with a mean wavenumber shift δ Δν av close to the 0.2 cm -1 overall shift (Figure 8) and a very constant standard deviation. This proves that the stress field produced by the in-pore pressure change is superimposed on an initial spatial stress/strain distribution observable in the Raman map under LV conditions (Figure 7, left). The southern zone shows a different behavior with a higher Raman shift difference (0.249 cm -1 ) corresponding to a Raman pressure of 28 MPa. It is worth noting that the highest shift difference represented by the red/orange pixels in Figure 9 follows the crystallographic axis.

Nonvisible Failure of the Container (CsCl Inclusion

). The experiment was carried out repeatedly (see Section 2.4 with the cyclical procedure). After around 80 cycles carried out over 3 weeks, with tension lasting from some minutes to some days during each cycle, a brutal change of T H was observed, expressing an internal failure (usually called "inclusion stretching" in the inclusionist's glossary; see ref. [START_REF] Bodnar | Systematics of stretching of fluid inclusions. I. Fluorite and sphalerite at 1 atmosphere confining pressure[END_REF] for a review about the different types of damage in inclusions) and a change of the cavity volume. T H was then upshifted by 16 °C, and logically, the T IN changed by a similar value, at T IN = 111 °C (Figure 6). Despite this thermometric evidence, the change in the cavity shape or volume was so small that it could not be observed optically (this is usual with this "inclusion stretching"), leading us to infer that the volume change resulted from the formation of damage at the microscale as a network around the inclusion, below the optical limit.

The Raman maps were then measured at 140 °C, that is to say, for the same T-T IN difference as for the previous data sets.

They showed a very weak Raman shift between the biphasic f10 and the monophasic situations (Figure 10), much weaker than in the precracked state (Figure 8). This fact matches very well with the presence of a dense network of microcracks, diminishing the lattice stress detectable with Raman spectroscopy.

ACS Earth and Space Chemistry

It means that, despite the low stress thus registered, the 549 cyclical procedure led to a real damage all around the cavity 550 materialized by a brutal increase of T H . The local damage is 551 also evidenced by the Raman maps recorded after the damage 552 of the inclusion, which show much lower values of the peak 553 stress (i.e., for the largest tension in the inclusion), +6.7 MPa 554 (Figure 10), compared to the precracked state (+20 MPa in 555 the plane, Figure 8). This large drop of the in-plane stress is 556 interpreted as the presence of a dense network of microcracks, 557 able to adjust the liquid tensile stress by a multitude of 558 microdisplacements inside the crack, instead of provoking a 559 lattice stress detectable with Raman spectroscopy. 560 3.3. Fracturing the Container (NaOH Inclusion). The 561 NaOH sample contains four inclusions shaped as "negative f11 562 quartz crystals" (Figure 11), numbered 1 to 4. T H is revolving 563 around 140 °C, corresponding to a mean water density of D ≈ 564 926 kg/m 3 according to IAPWS EoS (pure water approx-565 imation 37 ) (Figure 11). Inclusions 1-3 did not show any 566 instantaneous nucleation in the positive (Celsius) temperature 567 area. At RT and after isochoric cooling from T H , the liquid 568 tension in the three monophasic inclusions 1-3 was -132 569 MPa (Figure 11). Inclusion 4 had T IN equal to 44 °C so this 570 inclusion was never under tension under ambient conditions 571 while the three others were, giving us a "reference" inclusion 572 without a fluid pressure effect on the nearby surroundings.

573

With this RT superheating installed in inclusions 1-3, we 574 started observing the inclusions subjected to a permanent 575 tension from the in-pore liquid. Surprisingly, after only some 576 weeks and before any Raman records or accumulation of data 577 on the inclusions of interest, the host quartz cracked all around the three inclusions hosting tensile liquid, consistent with the crystallographic axis (Figure 11, bottom-right panel).

The quartz failure occurred with optically visible cracks appearing at the edge of the three inclusions 1-3. The bright border is underlined by a similar black edge when zooming: the cracks cannot be due to an optical artifact. Inclusion 4, which was never subjected to long tension (T IN > RT), was not damaged. The coexistence between the two mechanical behaviors, so well correlated with the microthermometric features, is a clear indication that the liquid tension is able to fracture the surrounding quartz. 12). In the equatorial plane, the compressive axial stress □ z is the largest contribution to the total pressure while the magnitudes of the tangent and radial stresses □ t and □ r are much smaller.

The tangent stress σ t is compressive as well and contributes to a lesser extent to the total pressure. On the other hand, radial stress σ r does not have a monotonous change. It has a negative value in the very vicinity of the edge, quickly inhibited by the global flattening of the inclusion yielding a compressive stress.

Real Inclusion Shape and the Stress Concentration

in the CsCl Inclusion. The Raman shift can be evaluated from either from Gruneisen's constant as shown above or from the deformation phonon potentials, usually noted p, q, and r,

ACS Earth and Space Chemistry

617 which are the components of the fourth rank tensor relating 618 the changes of the spring constants of the lattice and the 619 deformation tensor. [START_REF] Anastassakis | Effect of static uniaxial stress on the Raman spectrum of silicon[END_REF] They generally are determined from 620 hydrostatic pressure Raman experiments using diamond anvil 621 cells, uniaxial stress experiments, or strained thin film grown by 622 strained-layer epitaxy. It is worth noting that the third potential 623 r relates the shear strain components of the tensor. Briggs and 624 Ramdas 63 determined the deformation potential constants of 625 the zone-center optical phonons of α-quartz under uniaxial 626 stress experiments. Finally, the Raman shift can also be 627 evaluated from the stresses using the components of the fourth 628 rank elastic tensor relating stress and strain.

629

In the present case, the radial, tangent, and axial stresses are 630 also the principal stresses because of the symmetry of the 631 geometry and the loading, and it is therefore possible to derive 632 the expected Raman shift from either the principal stresses or 633 the radial, tangent, and axial stresses. Here, the constants 634 determined by Briggs and Ramdas [START_REF] Briggs | Piezospectroscopy of the Raman spectrum of -quartz[END_REF] were used to calculate the 635 464 cm -1 Raman shift from the stresses considering that the z-636 axis of the model coincides with the crystallographic c-axis as 637 shown in Section 2.1: 646 where P is in MPa and δν in cm -1 . 647 Both Raman shifts were evaluated along with the radial axis f13 648 (Figure 13). The very good agreement between Briggs and 649 Gruneisen formalisms confirms that the Gruneisen's constant-650 based approach was really reliable in this study, as proposed in 651 Section 2.3. 652 Focusing on the inclusion described in Section 3.1, the 653 Raman pressure at the bottom-right corner of the inclusion 654 (Figure 9) was compared with the simulation of the oblate f14 655 cavity with a shape factor b/a = 0.1 and anisotropic properties f14 (Figure 14). The simulated pressure (blue line FEA) quickly tends asymptotically to zero upon moving away from the inclusion, as expected from a direct application of the theory of solid mechanics. In the meantime, the experimental values showed that the sample is under stress over the total area (red stars, Figure 14; see also the green and blue pixels, Figure 9). The pressure derived from the Raman measures shows a slight decrease, especially close to the inclusion (gray area), where the values are in the same order of the simulation. Considering the difficulty in localizing the edge of the Raman map and the assumptions for the Raman pressure calculation and the simplified geometry for the simulation and despite the dispersion of the measured values, measurements and simulation can be estimated as consistent at least from a qualitative point of view. However, the Raman pressure seems to remain constant around 25 MPa when looking over the whole map even farther from the inclusion. This asymptotic mean pressure is not at all expectable and constitutes one of the main surprises of this experimental study. Although we do not have explanations for this bias, it may be compared with the Raman gap averaged over the whole map (0.212 cm -1 , i.e., 24 MPa). Another point to be considered is the presence of another cavity, which might not be alone, close to the studied inclusion. The stress fields induced by these surrounding cavities might produce such a general overstress. Finally, the general overstress may also be related to the localized strain/ stress induced by the inclusion synthesis procedure.

δν σ σ σ = • × + + • × - - - 0
The Raman pressure appears to be a good probe of the The pressure reaches its highest values rather in the middle of

ACS Earth and Space Chemistry

Focusing on the symmetry plane A (equatorial plane) that 722 contains edges B-C-D of the inclusion, we observe that two 723 principal stresses lie in plane A and that last one is 724 perpendicular to the plane. This can be explained by the 725 geometrical symmetry, quartz trigonal symmetry, and the 726 orientation of the inclusion with respect to the crystallographic 727 axes. The highest compressive stress is perpendicular to the 728 equatorial plane aligned with crystallographic c→ axis and is f16 729 very concentrated along the cavity edges (Figure 16). The 730 second compressive stress lies in the equatorial plane with 731 lower magnitude (Figure 16b, compare the color scale with 732 that of Figure 16a). The third principal stress also lies in the 733 equatorial plane and shows negative values corresponding to 734 tensile stress. It reaches its maximum magnitude around the 735 apex C (Figure 16c, tensile stress). The tensile stress direction 736 being along the line with the same origin and orientation as 737 experimentally observed in the studied inclusions (Figure 11) 738 confirms the possibility of crack formation along this line.

DISCUSSION

739 4.1. The Raman Pressure. Micro-Raman was here used to 740 dynamically follow the internal pressure in a crystal hosting a 741 variable pore tension. It is assumed that the applicability of the 742 Gruneisen formalism allows us to translate the Raman shift 743 into pressure. The anisotropy of monocrystalline quartz and the imperfections of the crystal lattice even under room conditions (see the Raman "irregularities" in Figure 7 left) limit this applicability. However, our study relies on the comparison of two datasets obtained at two different pore pressures, with all other parameters being constant. In addition, the shift is observed for two different Raman vibration modes, and it is very comparable. In that sense, the recorded Raman shift cannot be related to mechanisms other than lattice deformation and can be interpreted with confidence using the solid's linear elasticity.

Raman microspectrometry of silica minerals is commonly used to determine the pressure conditions of polymineralic assemblages containing various solid inclusions in minerals (e.g., refs [START_REF] Hemley | Pressure dependence of Raman spectra of SiO 2 polymorphs: -quartz, coesite and stishovite[END_REF][65][START_REF] Korsakov | Raman mapping of coesite inclusions in garnet from the Kokchetav Massif (Northern Kazakhstan)[END_REF][START_REF] Ashley | Geothermobarometric history of subduction recorded by quartz inclusions in garnet[END_REF][START_REF] Ashley | Significance of "stretched" mineral inclusions for reconstructing P-T exhumation history[END_REF][START_REF] Ashley | Quartz-in-garnet inclusion barometry under fire: Reducing uncertainty from model estimates[END_REF][START_REF] Cisneros | Evaluation and application of the quartz-inclusions in epidote mineral barometer[END_REF] ). The dominant compression mechanism is considered to be the reduction of (Si-O-Si) (e.g., ref. [START_REF] Parkinson | Present-day ultrahigh-pressure conditions of coesite inclusions in zircon and garnet: Evidence from laser Raman micro-spectroscopy[END_REF] ). A similar Raman barometry is used for the P-T conditions of the diamonds' source regions, based on determination of the internal pressure in olivine inclusions and the stresses in the surrounding diamond (e.g., refs. [START_REF] Izraeli | Raman barometry of diamond formation[END_REF][START_REF] Nasdala | Growth zoning and strain patterns inside diamond crystals as revealed by Raman maps[END_REF] ). More recently, the Raman shifts of sapphires around a zircon inclusion were converted to hydrostatic pressure and deviatoric components of the stress tensor, combining micro-Raman and photoluminescence imaging techniques, linked by the piezospectroscopic effects on Raman and photoluminescence spectra. [START_REF] Noguchi | Imaging of internal stress around a mineral inclusion in a sapphire crystal: application of micro-Raman and photoluminescence spectroscopy[END_REF] The internal stress was highly concentrated at the tips of the zircon crystal, resulting in a compressive stress of several hundreds of MPa and leading to most inclusions in radial ACS Earth and Space Chemistry 771 cracks propagating in the host sapphire crystal from the apex of 772 the hexagonal zircon. The present paper affords evidence that 773 cracking can be invisible while leading to significant variation 774 of the residual pressure. It could explain the growing 775 observations in natural samples of a wide scattering of residual 776 pressure sealed in inclusions (e.g., refs. [START_REF] Enami | Laser Raman microspectrometry of metamorphic quartz: A simple method for comparison of metamorphic pressures[END_REF][START_REF] Kouketsu | Evaluation of residual pressure in an inclusionhost system using negative frequency shift of quartz Raman spectra[END_REF][START_REF] Taguchi | Metamorphic record of the Asemi-gawa eclogite unit in the Sanbagawa belt, southwest Japan: Constraints from inclusions study in garnet porphyroblasts[END_REF][START_REF] Zhong | Zircon and quartz inclusions in garnet used for complementary Raman thermobarometry: application to the Holsno̷ y eclogite, Bergen Arcs, Western Norway[END_REF] ). 777 4.2. Crystal Host Stress Field. The comparison between 778 simulation and Raman measurements is to some extent 779 hindered by the heterogeneity displayed by the sample. 780 When recording the spectroscopic signatures at RT, namely, 781 the reference state without any external constraint, the data set 782 showed significant variations of the Raman peaks across the 783 whole map. Raman peaks showed a distribution into four 784 quadrants, which remained almost unchanged whichever the 785 fluid pressure (Figure 5). These measurements demonstrate 786 that the host crystal is not the perfect, stress-free monocrystal 787 assumed in the simulations. Experiments confirm the idea that 788 the elastic deformations induced by the growing fluid pore 789 tension retain and even magnify the initial imperfections. The 790 Raman maps clearly show that the tension-driven shifts are 791 amplified along the crystallographic axes. Therefore, both the 792 Gruneisen formalism to convert Raman signatures into 793 pressure and a direct comparison of such values to FEA 794 simulations done through isotropic elastic theory can only 795 approximately describe the complex mechanical behavior of 796 quartz.

797

In spite of these limitations, the Raman measurements do 798 retrieve features expected from elastic theory: (i) a 799 compressive stress develops within the equatorial plane and 800 (ii) the stress concentration is maximal in the vicinity of the 801 cavity tip (Figure 14). Furthermore, the shape of the pore 802 cavities plays a role in the stress field to build up a potentially complex 3D result (Figures 15 and16). This complexity is certainly increased by the pre-existence, before installing any new constraints, of a certain degree of heterogeneity that is visible in the Raman maps at RT (Figure 7 [START_REF] Vityk | Fluid inclusions as tectonothermobarometers: Relation between pressure-temperature history and reequilibration morphology during crustal thickening[END_REF] ). In a later study, cracking was observed when internal underpressure only exceeded 270 MPa. [START_REF] Vityk | Textural evolution of synthetic fluid inclusions in quartz during reequilibration, with applications to tectonic reconstruction[END_REF] In our case, the amplitude of tensile stress in the inclusion, analogous to an underpressure stress, is either ∼200 MPa

ACS Earth and Space Chemistry

(CsCl sample) or ∼140 MPa (NaOH sample). In spite of internal underpressure lower than [82-83]'s, we observed the formation of cracks with the lower differential stress. The simplest explanation for the difference is the effect of confining pressure, absent in our present study, while it is known to inhibit cracking. [START_REF] Vityk | Textural evolution of synthetic fluid inclusions in quartz during reequilibration, with applications to tectonic reconstruction[END_REF] An additional effect lies in the geometrical shape of the fluid inclusion, commonly spheroidal for natural inclusions. Our FEA modeling shows that the stress concentration at the tip of the inclusion reached a larger amplitude for the negative crystal shape (as displayed by the SFIs we used) than the spheroidal one. The geometry thus promotes the concentration of stress in the very limited area of the apex, which may facilitate crack initiation. Furthermore, we can speculate a little more about the micromechanical mechanisms that underlie the development of cracks at low levels of tension. For instance, the microprestress theory 4 offers a comprehensive view of how capillary and disjoining pressures can drive the nucleation and growth of tensile microcracks in the complex pore network of an aging cement. In particular, hindered adsorption inside the nanoscale pores (1 to 10 molecular width and 0.3 to 3 nm) within hardened cement paste plays a pivotal role. It is appealing to hypothesize, for our experiments, the existence of tiny quantities of water present along the quartz dislocations comparable to the gel pores for cement. The liquid tension existing in the cavity might therefore spread over the nearby crystal lattice owing to minute water-bearing dislocations, propagating the tensile stress farther from the cavity and concentrating the constraints.

CONCLUSION AND PERSPECTIVES

The tension of water occluded in a crystalline cavity has a pore-scale micromechanical effect, as exemplified by the visible cracking for the NaOH sample, and the microthermometric features change (T H increasing) for the CsCl sample. We demonstrated here that Raman spectroscopy is able to monitor in space and time the pressure in the quartz crystal around a fluid inclusion. The Raman maps revealed that the stress fields induced by the pore tension within the crystalline host rock are nonuniform and concentrated along the crystallographic orientations and, probably, crystal defects. The exact value of stress at these overstressed sites cannot be ascertained by the Raman technique because they may operate at a much smaller scale than the Raman scattering volume. The numerical-model simulated stress field and stress concentration match the measured 3D distribution of stress well and outline the role of geometry in terms of stress concentration in the region where fracturing was observed. On the contrary, the Raman estimates of how stress changes with distance from the inclusion do not match the elastic theory, showing that the crystal defects probably play a role in the stress distribution. Technically, micro-Raman associated with SFI microthermometry appears to be a powerful tool for deciphering the micromechanical mechanisms involved in complex field stress in rocks and porous media.

To summarize, our study demonstrates unambiguously that a permanent low stress as well as repeated cyclic stresses established in a crystal lattice by reducing the in-pore pressure in an occluded water-bearing cavity (related to capillary-type effects) is able to damage the hosting quartz despite its high resistance. The mechanical role pressurized pore fluids can play in weakening and embrittling rocks has long been known (e.g., refs. [START_REF] Paterson | Experimental rock deformation -The brittle field[END_REF][START_REF] Bruno | Pore pressure influence on tensile fracture propagation in sedimentary rock[END_REF] ), and this paper is a first step in envisioning the role of

Figure 1 .

 1 Figure 1. Photomicrographs and EBSD measurements of the two inclusions of interest (L = liquid aqueous solution; V = vapor). A. CsCl inclusion entirely filled by aqueous solution. B. NaOH inclusion. For each inclusion, two stereoplots show the orientation of the [0001] and [112̅ 0] axes within the host quartz, located perpendicular and parallel to the polishing plane, respectively. The redundant coordinate system (a 1 , a 2 , a 3 , c) is used (hexagonal lattice).

  software (manufacturer EDAX, Mahwah, USA), available at 161 BRGM (Orleáns, France). The working distance was of about 162 18 mm, at an accelerating voltage of 25 kV. 163 The samples are platelets containing the a-and b-axes and 164 thus perpendicular to the crystallographic c-axis (Figure 1). 165 2.2. Method 1: Establishing Tensile Conditions in the 166 Pore Fluid. High tension can be produced in a liquid 167 occluded inside a closed cavity, by making it monophasic 168 through increasing temperature and then cooling down this 169 monophasic liquid below the saturation curve along a constant 170 volume trajectory. This procedure was largely used in the past 171 to study the phase diagram of tensile water (e.g., refs. 27-31 ) 172 and the corresponding changes of phase equilibria (e.g., 173 refs. 33,34 using fluid inclusions; refs. 35,36 using Berthelot tubes). 174 The starting pore material generally contains a biphasic 175 liquid + vapor assemblage (L + V). With progressive heating, 176 the liquid density and bubble pressure change, and the trapped 177 liquid progressively invades the whole inclusion space (A to B 178 f2 path, Figure 2) up to a particular density where there is no 179 more vapor (point B, Figure 2). The inclusion is then 180 recognized to be homogenized at the corresponding 181 "homogenization temperature" (T H ). Further heating drives 182 the P-T conditions in the inclusion along the isochore curve 183 inside the stable domain of liquid (B to C path, Figure 2). The 184 second step is the progressive cooling of the sample, which 185 follows the isochoric path as long as the inclusion remains 186 homogeneously filled with liquid (C to D path without 187 nucleation at B, Figure 2). Once T H is crossed without bubble 188 nucleation, the monophasic liquid becomes metastable; that is, 189 it is less stable than the liquid-vapor assemblage. The 190 nucleation does not occur as long as the energy cost of the 191 vapor-liquid interface fabrication is higher than the gain 192 afforded by the bulk energy of the most stable phase. 31,38 Such 193 a monophasic liquid is superheated (located to the right of the 194 saturation curve) and follows an isochoric cooling down to a 195 certain temperature at which a bubble appears suddenly. This 196 is the instantaneous nucleation temperature (T IN ), which is 197 always located within the tensile domain: nucleation in fluid 198 inclusions disobeys the saturation conditions. After nucleation, 199 the newly formed LV assemblage moves instantaneously back 200to the saturation curve conditions (D to E path, Figure2). This

Figure 2 .

 2 Figure 2. Phase diagram of pure water in the P-T space, calculated with the IAPWS-95 equation of state.[START_REF] Wagner | The IAPWS Formulation 1995 for the Thermodynamic Properties of Ordinary Water Substance for General and Scientific Use[END_REF] Each isochore is labeled with the molar volume in cm 3 •mol -1 . The CsCl inclusion micrographs illustrate the fluid content present in the cavity at each characteristic step.

  357 factor b/a. For a sphere, b/a = 1, the axial and tangent stresses 358 are equal and compensate the radial stress, and the resulting 359 pressure is null. Each stress does not exceed the pressure, and 360 no stress concentration effect occurs. On the contrary, when 361 the spheroid is flattened, the axial and tangent stresses are not 362 equal any longer, and they do not compensate the radial stress. 363 The axial stress largely predominates the other stresses and 364 reaches up to 12 times the fluid inclusion pressure. The 365 positive values of the axial stress and pressure indicate a 366 compressive state in equatorial plane although a tensile stress is 367 f4

368

  This effect is due to the flattening of the spheroid and is 369 magnified when the shape factor decreases. 370 Once validated, the ellipsoid-based model was used to 371 simulate the experiments. A negative (tensile) pressure was 372 applied on the 2D elliptical edge to simulate the in-pore 373 pressure change. No gradient or interfacial effects of any sort 374 that may "buffer" the liquid-to-solid pressure transmission were 375 taken into account. The anisotropy of α-quartz was considered, 376 and the elastic moduli were defined according to Calderon et 377 al. 57 378 Second, the sample was modeled with a 3D prismatic solid 379 containing a faceted inclusion to look like the inclusion 380 f5 observed with Raman spectroscopy (Figure 5). Only an eighth 381 of the sample was modeled considering symmetry conditions, 382 and the dimensions were defined to match the inclusion 383 actually considered. Zero displacement was applied to the 384 three symmetry planes (SP), and a negative pressure was again 385 applied on the inner surfaces of the inclusion. The same

Figure 3 .

 3 Figure 3. Normalized stress along the radial axis. The normalized stresses are calculated by dividing the stresses by the pressure applied onto the inclusion wall. The normalized radial coordinate is the radial coordinate divided by the half-length of the spheroid, r* = r/a. The stresses denoted "Esh" were calculated following Eshelby's theory 45 by applying script developed by Meng et al. 51-53 and according to the superposition principle for linear elastic mechanics (top-left of the figure).Stresses that were noted by FEA were calculated using Comsol Multiphysics software. For all the calculations, Young's modulus and Poisson's coefficient were 100 GPa and 0.06, respectively, according to ref.56

3883. 1 .

 1 Stressing the Container (CsCl Inclusion). The CsCl 389 fluid inclusion has a T H of 205.6 ± 0.4 °C and T IN of 89 ± 1.4 390 °C. Drawing a straight line from the synthesis conditions (630 391 MPa and 480 °C) to the (T H ,P H ) pair (dP/dT = 2.3 MPa/°C), 392 and extrapolating this line to the T IN value, gives the liquid 393 tension existing in the inclusion just before the nucleation, at f6 394 -266 ± 3 MPa (Figure 6).

  395

  spots correspond to regions where the Raman line intensity decreases below a threshold and were used to locate the inclusions. The barycenter of the black spot corresponding to the inclusion of interest was used as the origin of the map. Finally, a micrograph of the inclusion of interest has been superimposed by aligning the lattice directions determined with EBSD with the Raman map with a perfect match. The dark spot is slightly bigger than the inclusion image because of the empirical determination of the threshold. Applying the Schmidt and Ziemann 42 relation at 120°C gives an expected wavenumber of 462.4 cm -1 , which can favorably be compared with our mean measured wavenumber 462.8 cm -1 . Because of uncertainties (absolute calibration of two different spectrometers, accuracy on the wavenumber of a "reference" quartz totally stress-free at RT), these values look similar. Though the wavenumbers of the map seem to follow a rather narrow spread, their spatial distribution is well organized around the inclusion. Four zones aligned with the crystal axes show a remarkable symmetry (Figure7, left). In the upper and lower triangular zones (where green predominates), the wavenumber is close to 462.75 cm -1 while it is close to 463 cm -1 in the lateral zones (orange/red pixels). The sharp variation in internal pressure between these domains leads us to think of some localized strains on planar defects, presumably connected with how the inclusion formed (thermal fracturing followed by SiO 2 overgrowth).The map obtained from inclusion in the monophasic state for the same A 1 Raman mode is presented in Figure7on the right. Again, the micrograph of the inclusion of interest has been superimposed. In this figure, the dark spot is hidden by the inclusion image, indicating that the dark spot is smaller than that of the biphasic state map. That is also the case for the other spot located at the bottom-right of each map. This shows that the confocal z plane slightly changed. The mean wavenumber is 463.0 cm -1 , ranging from 462.8 to 463.3 cm -1 . Again, the temperature-corrected wavenumber calculated with the Schmidt and Ziemann 42 relationship compares favorably with the "monophasic" average value. However, it is the gap between the "biphasic" and "monophasic" values which is of special interest. It is very small (0.25 cm -1 between 462.75 and 463 cm -1 ), but both measurements were performed on the same spectrometer with the same calibration procedure, and on the same points of one quartz sample. Only the spectrometer stability is involved here. It is worth noting that the color scale for the monophasic state map slightly differs from the LV state one. The color scale of each map has been chosen to emphasize the similarities between the spatial distributions of the colors and consequently of the wavenumbers. Globally, the increase of the Raman wavenumbers of the monophasic state map indicates the development of compressive stress although no optical evidence of quartz failure was recorded during this experiment, which therefore follows a purely elastic behavior.

Figure 4 .

 4 Figure4. Simulation of a host with a spherical or spheroidal inclusion with a negative (tensile) pressure applied on the wall cavity. See Figure1for more details on the model. Normalized pressure simulated with FEA and depending on shape factor b/a, maps, and values along the radial axis. For each shape factor, the continuous line and marks represent the FEA-based pressure and the Eshelby-based pressure, respectively.

Figure 5 .

 5 Figure 5. Geometry of the 3D mechanical model of the inclusion for the FEA. Because of symmetries (SP), only an eighth of the solid is modeled. The image size is not scaled with the prism size for visibility purposes.

Figure 6 .

 6 Figure6. Thermal procedure, microthermometric features, and schematics of the cyclic tension applied to the sample, before (label 1, Section 3.1) and after (label 2, Section 3.2) the stretching irreversible failure of the inclusion.

Figure 7 .

 7 Figure 7. Raman map of the wavenumber of the high wavenumber A 1 quartz mode, with infilling fluid at the liquid-vapor saturation state (left) and monophasic liquid state (right). The wavenumber values are given by the color scale. The inclusion of interest has been superimposed as well as the crystal axes determined with EBSD.

Figure 8 .

 8 Figure 8. Raman line center histograms at 120 °C of the 464 (left) and 206 (right) cm -1 bands of quartz around a fluid inclusion filled with either a biphasic liquid-vapor stable mixture (green, blue = repetition to test for reproducibility) or a monophasic tensile liquid (red).

Figure 9 .

 9 Figure 9. Shift difference between biphasic and monophasic modes for the 464 cm -1 Raman mode. Top-left: Map of the shift difference. Top-right: global histogram of this map. Bottom: Histograms of the four regions indicated by the red rectangles on the map, with a geographic-like denomination.

Figure 10 .

 10 Figure 10. Histograms of the line center of the high wavenumber band, at 140 °C, along the equatorial plane of the inclusion, after the T H -recorded damage.

3. 4 .

 4 Numerical Models of Stress Concentration and Fracturing. 3.4.1. Ellipsoidal Model and Pressure in the Equatorial Plane. The simulated stresses and pressure resulting from FEA simulation based on oblate spheroids f12 with a shape factor b/a = 0.1 are represented in Figure 12. The pressure is highly concentrated close to the inclusion edge where it reaches four times the fluid inclusion pressure and decreases sharply along the radial axis. The pressure distribution corroborates the Raman observations (Figure 10): a compressive stress (positive) lies within the equatorial plane around the cavity (red color). This can be explained by the flattening of the spheroid, like what occurred with the remote stress above. The simulation indicates that the principal stresses are aligned with the geometrical axes, that is, the radial, tangent, and axial directions of the ellipsoid (see the r, t, and z directions, respectively, in the inset top-right in Figure

Figure 11 .

 11 Figure 11. Upper panel, installing the "enduring tension" situation by isochoric cooling of a high-density SFI. Bottom left, observation of three fluid inclusions numbered (1-3) with microcracking in the quartz at each corner of the negatively shaped cavity. The fourth (num. 4), which does not tolerate tensile stress (see the bubble), is intact. Bottom-right, crystallographic axes measured by EBSD are reported on the right micrograph.

Figure 12 .

 12 Figure 12. 2D axisymmetric mechanics model of an ellipsoidal-shaped inclusion. Bottom: Normalized pressure map around an inclusion. The vertical symmetry axis is the revolution axis. Top: Stresses and pressure magnitude along the radial axis in the optical plane.

  stresses are in MPa and δν in cm -1 . The Briggs' 639 potentials were determined with an uncertainty of ±0.02 • 10 -2 640 cm -1 /MPa. Note that they are here positive because of the 641 change in the sign of the stresses. The uncertainty of the 642 Raman shift evaluated from uncertainty propagation depends 643 on the stresses and reaches a maximum value of 0.5 cm -1 .644The Raman shift was also calculated following the Gruneisen 645 formalism:

  stress fields around the water-bearing cavity, and a highly resolved measurement of the instantaneous/local stress is certainly possible according to our preliminary results. It is also clear that any change in the in-pore pressure drives a significant modification in the poromechanical stress field, over an unexpected long-range scope. To address these findings, further studies should involve the anisotropic strain propagation in the matrix along preferential axes (crystallography, dislocations, etc.) by carefully plotting how the quartz band varies from place to place, at saturation pressure. It was already clear (Figure 7 left) that our hydrothermal monocrystalline quartz was not uniform in terms of Raman shift. 3.4.3. Real Inclusion Shape and the Fracturing of the NaOH Inclusion. The objective of the FEA simulation performed on the 3D faceted cavity (Figure 5) was to model the fracturing of the NaOH hexagonal shaped inclusions as a result of tensile liquid in the cavity (Figure 11). As a preamble, it is well known that geometrical singularities yield artificially increased stress in finite elements depending on the meshing. In addition, the real interface between the inclusion and the contained liquid is not as clear as in the geometrical model where the interface corresponds to a wall. Therefore, values of the simulated pressure must be carefully considered. The main purpose is to obtain tendencies that give us a better understanding of crack formation. As expected, the (negative) tensile pressure applied on the inner wall of the cavity (-200 MPa) yields a tensile state within the major part of the solid domain proved by the negative value of the simulated pressure f15 (Figure 15, left). Like for the oblate spheroid model, the shape factor effect is observed as the pressure increase concentrates around the inclusion, reaching positive values along the cavity edges lying in the symmetry plane A, significantly higher than the applied inner pressure and indicating a compressive state in the plane induced by the predominant flattening of the cavity.

Figure 13 .

 13 Figure 13. Raman shift evaluation along the radial axis.

Figure 14 .

 14 Figure 14. Left, Zooming in the southern area having the larger average wavenumber Raman shift (i.e., Figure9). Right, Raman pressure (red stars) vs simulation (blue curve). The behavior close to the interface is highlighted as a gray area.

Figure 15 .

 15 Figure 15. Pressure field showing (left) stress concentration along the longest edge lying in the symmetry plane XY (see Figure 5) near inclusion wall angles and (right) domains of negative/positive stress in the XZ symmetry plane.

  left). 4.3. Stress Accumulation and Fracturing Threshold. It is difficult to be quantitative here; however, the stress generated around the water-bearing cavity appears moderate, especially compared to the observations of the damage (T H variations, apex fractures) in the host quartz. One proposition for reconciling low stress and evidence of cracking is the role played by fatigue effects in the resistance of materials. Stresses that vary with time and in intensity are able to cause the rock failure while below the yield strength, and even in the linear elastic domain (e.g., refs. 5,79-81 ). The best illustration is what is known as "reequilibration" in fluid inclusions studies, which is a modification of the inclusion shape and/or a fluid leakage because of the effective stress, namely, the pressure difference between host mineral stress and cavity-filling fluid pressure (e.g., refs 82-86 ). One strong hypothesis, often claimed to infer various deep fluid conditions from fluid inclusions, is their isochoric behavior along the rock pathway throughout the Earth's crust (burial and/or uplift). To validate it, the resistance of quartz to differential stresses was measured on SFIs for internal under-or overpressure (respectively P solid > P fluid or P solid < P fluid ). Overpressure causes reequilibration from reaching 100 MPa. 82,85 Underpressure caused visible cracks from around 200 MPa, while reequilibration proceeded by dissolution-precipitation in the range of 100-200 MPa (e.g., ref.

Figure 16 .

 16 Figure 16. Principal stress magnitude map near the corner of the symmetry plane XY (see Figure 5, plane noted A). (a) Principal compressive stress (σ ⊥ ) perpendicular to the plane. (b) Second principal compressive stress (σ II ) field near the corner of the symmetry plane XY (Figure 5, plane noted A). The double black arrows indicate the direction of σ II along the line corresponding to the potential crack direction. (c) Principal tensile stress (σ T ) field near the corner of the symmetry plane XY (see Figure 9, plane noted A). The double black arrows indicate the direction of σ I along the line corresponding to the potential crack direction. Note that negative/positive values indicate tensile/compressive stress.

  896 tensile pore fluids. Furthermore, this result poses the question 897 of the mid/long-term safety of geological storage (CO 2 , 898 nuclear wastes, energy) wherever the geological formations are 899 subjected to drying operations that make it possible to install 900 liquid tension in large masses of unsaturated hydrosystems. In 901 particular, the use of clayey formations to guarantee the 902 confinement of the nuclear species owing to this impermeable 903 barrier must be questioned if the consequent tension of tunnel 904 ventilation was able to enhance permeability in the massif by 905 creating fissures. In the same way, the integrity of the caprocks 906 in CO 2 storage is a crucial property that prevents any leakage 907 from occurring. However, the injection of supercritical CO 2 908 severely dries out the formation, especially at the highest part 909 where the CO 2 bubble accumulates at the very contact with the 910 caprock, establishing the right conditions for tension-driven 911 processes. 912 During the course of the study, the initiation of micro-913 cracking and brittle failure was also discussed in association 914 with the water distribution in rock. This aspect is of great 915 significance in natural systems, beyond the question of storage 916 safety. At high pressure and temperature, where deformation 917 proceeds mainly as a result of temperature-dependent viscous 918 creep processes (dislocation-or diffusion creep), fracturing is 919 still present. These brittle precursors promote grain-size 920 reduction and therefore the fluid flow and chemical reactions 921 that activate water corrosion. 89-91 High fluid pressures, up to 922 confining pressure values, are necessary to open cracks, but the 923 mechanisms are not well understood, even if we know that the 924 physical state of the water present in mineral assemblages with 925 very low porosity is often involved as a triggering 926 mechanism.[START_REF] Brok | Effect of microcracking on pressure-solution strain rate: The Gratz grain-boundary model[END_REF][START_REF] Nakashima | Infrared microspectroscopy analysis of water distribution in deformed and metamorphosed rocks[END_REF] This study demonstrates that reducing fluid 927 pressure also induces cracking mechanisms rooted at the small 928 scale, even if the generality versus specificity of such 929 mechanisms in the various pressure contexts is still to be -ISTO, UMR 7327, Univ. Orléans, CNRS, 934 BRGM, Orléans F-45071, France; orcid.org/0000-0002-935 3825-663X; Email: lionel.mercury@univ-orleans.fr

https://dx.doi.org/10.1021/acsearthspacechem.0c00224ACS Earth Space Chem. XXXX, XXX, XXX-XXX

https://dx.doi.org/10.1021/acsearthspacechem.0c00224ACS Earth Space Chem. XXXX, XXX, XXX-XXX D

https://dx.doi.org/10.1021/acsearthspacechem.0c00224ACS Earth Space Chem. XXXX, XXX, XXX-XXX F

https://dx.doi.org/10.1021/acsearthspacechem.0c00224ACS Earth Space Chem. XXXX, XXX, XXX-XXX L

Authors

937

Emmanuel De Bilbao -CEMHTI, UPR3079 CNRS, Univ.


Qishuai Liu
email: qishuai.liu@huskers.unl.edu

Emmanuel Moulay
email: emmanuel.moulay@univ-poitiers.fr

Patrick Coirault
email: patrick.coirault@univ-poitiers.fr

Qing Hui
email: qing.hui@unl.edu

Deep Learning Based Formation Control for the Multi-Agent Coordination

HAL is a

I. INTRODUCTION

The cooperative control of multi-agent systems (MASs) has drawn a great attention in recent years in different areas [START_REF] Ge | A survey on recent advances in distributed sampled-data cooperative control of multi-agent systems[END_REF] including reconnaissance, surveillance, and security [START_REF] Fields | Effective robot team control methodologies for battlefield applications[END_REF], [START_REF] Ding | Hybrid systems in robotics[END_REF]. In such applications, maintaining the network topology and connectivity of agents are important for some tasks including target localization, oceanic search, and undersea oil pipeline maintenance [START_REF] Rezaee | Pursuit formation of double-integrator dynamics using consensus control approach[END_REF]- [START_REF] Li | Receding horizon formation tracking control of constrained underactuated autonomous underwater vehicles[END_REF].

Formation control is one of the most interesting research topics among cooperative control of MASs [START_REF] Chen | Formation control: a review and a new consideration[END_REF]. In many applications, a group of autonomous robots are required to follow the specified trajectory and maintain the network topology at the same time because of the system cost, redundancy, and structure flexibility of the system [START_REF] Stilwell | Platoons of underwater vehicles[END_REF]- [START_REF] Robertson | Formation sensing and control technologies for a separated spacecraft interferometer[END_REF]. In formation control, for a group of coordinated agents, different control topologies can be adopted according to different tasks. Some of them may need a leader for other agents so that they can follow the leader in a specified way. Those robots carry onboard sensing equipment to maintain the relative topologies between them. Generally speaking, the robots have the limited communication ability, thus they cannot use all the global information of the other robots. Given this situation, the design of the controller for each robot should not be centralized but be based on the local information.

In order to design the distributed controller for the multi-agent formation, we need to consider some issues such as the stability of the controller, controllability of different formation patterns, and safety and uncertainties in the formation control. The traditional control approaches to solve these issues are to develop the leader-follower strategy [START_REF] Cowan | Vision-based follow-the-leader[END_REF], virtual structure approach, and behavior-based method. The formation control is classified as regulation control and tracking control in some situations. Generally, this may need to consider the system dynamic and design the controller carefully. In some cases, this is hard to implementation. In order to overcome this drawback, some researchers develop artificial neural network based controllers [START_REF] Cai | Intelligent Control: Principles, Techniques and Applications[END_REF]. Deep learning (DL), which is based on the nonlinear continuous functions [START_REF] Cybenko | Approximation by superpositions of a sigmoidal function[END_REF], [START_REF] Hornik | Multilayer feedforward networks are universal approximators[END_REF], is part of machine learning [START_REF] Ao | Machine Learning and Systems Engineering[END_REF] and is dedicated to the modeling of nonlinear behaviors by a learning process involving deep neural networks (DNN) [START_REF] Lecun | Deep learning[END_REF]. In the past several years, the deep learning method is applied in many aspects such as pattern recognition [START_REF] He | Deep residual learning for image recognition[END_REF], speech recognition [START_REF] Hinton | Deep neural networks for acoustic modeling in speech recognition: The shared views of four research groups[END_REF], and computer vision [START_REF] Janai | Computer vision for autonomous vehicles: Problems, datasets and state-of-the-art[END_REF].

One of the machine learning based methods is the adaptive dynamic programming, which has been used as the model-based algorithm to design the controller for dynamical systems [START_REF] Liu | Adaptive Dynamic Programming with Applications in Optimal Control[END_REF]. This method uses a cost function to find the feedback controller named policy by solving the Bellman equation of the value function [START_REF] Liu | Adaptive Dynamic Programming with Applications in Optimal Control[END_REF]. DNN approximations of value function and policy are used in ADP to stabilize discrete-time systems. The policy can be obtained from the value function [START_REF] Jiang | Data-driven adaptive dynamic programming schemes for non-zero-sum games of unknown discrete-time nonlinear systems[END_REF].

Another machine learning based method is reinforcement learning, which is used to develop the model-free algorithms of controllers for dynamical systems [START_REF] Fourati | Stabilization of unknown Fig. 7. The initial position for all robots. Fig. 8. The robots are running. Fig. 9. The final position for all robots. nonlinear systems using neural networks[END_REF]. The main aim of this method is to learn from experiments on what to do in different situations so that it can maximize a reward function [START_REF] Panait | Cooperative multi-agent learning: The state of the art[END_REF]. By using this method, the designed controller is a model-free strategy that can be obtained by solving the Bellman equation of the value function. This method can be used for both stochastic processes using MDP and deterministic processes. By incorporating the DNN into the reinforcement learning, 978-1-7281-0084-5/19/$31.00 c 2018 IEEE a new method named deep reinforcement learning (DRL) has been proposed for the control of dynamical systems [START_REF] Duan | Benchmarking deep reinforcement learning for continuous control[END_REF], [START_REF] Hafner | Reinforcement learning in feedback control[END_REF]. The value function and the policy are approximated by DNN in the DRL approach.

The main drawback of these two methods is that the Bellman equation is hard to solve. Thus, the neural network can be used as a controller to ensure the control of dynamical systems [START_REF] Psaltis | A multilayered neural network controller[END_REF]. In our paper, we extend the neural network to use the DNN as the controller for discrete-time dynamical systems and use it for the formation control. This method, first developed in [START_REF] Li | Neural network control of unknown nonlinear systems[END_REF] for shallow neural networks and in [START_REF] Fourati | Stabilization of unknown nonlinear systems using neural networks[END_REF] for DNN under a restrictive assumption on the system, is based on the ability of a DNN for self-tuning its weights by using a backpropagation algorithm. Thus, it can avoid the solving the Bellman equation in ADP or DRL. Also, this method can be used if the vector field of the discrete-time nonlinear system is differentiable. Thus if the discrete-time nonlinear system is unknown, it is still possible to use numerical differentiation methods for backpropagation.

This paper is organized as follows: Section II gives some preliminaries and Section III introduces the deep learning based controller design. In Section IV, we introduce how to transform the formation control problem into the feedback control so that we can use the proposed algorithm to solve it. We also have done the experiment to verify the proposed algorithm in Section V. Finally, we give a conclusion in Section VI.

II. PRELIMINARY

The notation used in here is fairly standard. Specifically, R denotes the set of real numbers, N denotes the set of nonnegative integers, R n denotes the set of n-dimensional real column vectors, R n×m denotes the set of n-by-m matrices, (•) T denotes transpose, and (•) -1 denotes inverse, respectively.

• denotes the Euclidean norm. diag(x) denotes a square diagonal matrix with the elements of vector x on the main diagonal. Let G = (V, E) denote a dynamic directed graph with the set of vertices V = {v 1 , v 2 , . . . , v N } and E ⊂ V × V representing the set of edges. A graph with the property that

(v i , v j ) ∈ E implies (v j , v i ) ∈ E is said to be undirected. The adjacent matrix A ∈ R N ×N associated with the directed graph G is defined by nonnegative adjacency elements a ij as A(t) = [a ij (t)]. We assume (v i , v j ) ∈ E(t) if and only if a ij = 1, (v i , v j ) ∈ E(t) if and only if a ij = 0, and a ii = 1 for all i ∈ N , where N = {1, 2, . . . , N } denotes the node index of G(t). The set of neighbors of the node v i is denoted by N i (t) = {v j ∈ V : (v i , v j ) ∈ E(t), j = 1, 2, . . . , |N |, j = i}, where |N | denotes the cardinality of N . The degree matrix of a dynamic graph G(t) is defined by δ(t) = [d ij (t)], where i, j ∈ {1, 2, . . . , |N |} and d ij (t) = |N | j=1 a ij (t), if i = j, 0, if i = j.
Moreover, we can define Laplacian matrix of the dynamic graph G as L(t) = δ(t) -A(t). If the graph is called strongly connected if there exists a path from any node to any other node in the dynamic graph.

Next, let K ⊂ R n be a compact set such that 0 ∈ K and called the learning state space. Consider the following system:

x(k + 1) = f (x(k), u(k)), k ∈ N, x(0) = x 0 ∈ R n (1) where x(k) = (x 1 (k), . . . , x n (k)) ∈ K is the state vector, u(k) = (u 1 (k), . . . , u m (k)) ∈ R m is the control input with m ≤ n, and f : R n × R m → R n a continuous function such that f (0, 0) = 0 and f = (f 1 , . . . , f n ) with f i : R n ×R m → R.
The objective is to ensure that the solutions of system (1) converge to the origin by using a DNN as a controller. To this end, we consider the extended system

x(k -s + 1) = f (x(k -s), u(k -s)), k ∈ N, 0 ≤ s ≤ r (2)
where (x(k -r + 1), . . . , x(k + 1)) ∈ K r+1 is the state vector with

x(k -s + 1) = (x 1 (k -s + 1), . . . , x n (k -s + 1)) ∈ K, (u(k -r), . . . , u(k)) ∈ R m(r+1) the control input with m ≤ n and u(k -s) = (u 1 (k -s), . . . , u m (k -s)) ∈ R m ,
and r ≥ 0 such that the following assumption holds:

Assumption 1. Let E(k) = 1 2 x(k + 1) 2 (3) if lim k→∞ ( ∂E(k) ∂u1(k-r) , . . . , ∂E(k) ∂um(k-r) , . . . , ∂E(k) ∂u1(k) , . . . , ∂E(k) ∂um(k) ) = 0 then we have lim k→∞ x(k) = 0.
Next, let us recall the definition of a trajectory of the system (1) given in [29, p. 13]. Definition 1. Given an initial value x 0 ∈ K and a control sequence u(k), we define a trajectory x u (k, x 0 ) of the system (1) by x u (0, x 0 ) = x 0 and

x u (k + 1, x 0 ) = f (x u (k, x 0 ), u(k)), k ∈ N
The classical feedback control aims at finding u(x(k)) to stabilize the system (1), as shown by the Lyapunov based methods in [START_REF] Lin | Feedback stabilization of general nonlinear control systems: A passive system approach[END_REF], the model predictive control in [START_REF] Allgöwer | Nonlinear Model Predictive Control[END_REF] or the neural-network-based ADP algorithm in [START_REF] Sarangapani | Neural Network Control of Nonlinear Discrete-Time Systems[END_REF], following the framework illustrated in Figure 1.

III. METHOD

In this section, we will introduce the backpropagation algorithms used in the DNN. The DL feedback control aims at finding a feedback control U : K → R m as a DNN to ensure the convergence of the state of the system (1) toward zero by using the backpropagation algorithm and without using classical feedback control theory or the Bellman equation.

Let us denote U (k -s) = (U 1 (k -s), . . . , U m (k -s)) ∈ R m the output of the DNN with the input x(k -s) ∈ K where 0 ≤ s ≤ r. It means that the DNN as the r + 1 inputs x(k -s) ∈ R n and the r + 1 outputs U (k -s) ∈ R m where 0 ≤ s ≤ r. We are facing problems for backpropagation due to the fact that we have no target for the output U (k-s) of the DL feedback control and moreover U (k -s) ∈ R m whereas x(k + 1 -s) ∈ R n with m ≤ n. The main goal is to make x(k) converge to zero, so we consider the backpropagation error (3) for the input x(k -s) of the DNN. The tuning of the weights w ij (k) between neurons i and j for the input x(k -s) of the DNN is done by using the backpropagation algorithm

w ij (k -s + 1) = w ij (k -s) -λδw ij (k -s) = w ij (k -s) -λ ∂E(k) ∂w ij (k -s) (4)
with λ > 0 being the learning rate and the chain rule:

∂E(k) ∂w ij (k -s) = ∂E(k) ∂σ j (k -s) ∂σ j (k -s) ∂z j (k -s) ∂z j (k -s) ∂w ij (k -s)
where k ∈ N, 0 ≤ s ≤ r, and

σ j (k -s) = ϕ j (z j (k -s)) = ϕ j l w lj σ l (k -s)
in which k ∈ N and 0 ≤ s ≤ r. The only difference with usual backpropagation is the first term ∂E(k) ∂Uj (k-s) due to the fact that the backpropagation error E(k) does not depend directly on U j (k -s). Consider the following assumption for f :

Assumption 2. The function (x, u) → f (x, u) is differentiable on K × R m .
By using the system (2) and Assumption 2, we can obtain the following derivative for the backpropagation error (3) with respect to the output

U j (k -s), k ∈ N and 0 ≤ s ≤ r ∂E(k) ∂U j (k -s) = n i=1 x i (k + 1) ∂f i (x(k), U (k)) ∂U j (k -s) . (5)
In order for the backpropagation algorithm (3)-( 5) to converge, we must add the following assumption:

Assumption 3. The function (x, u) → f i (x, u) 2 is convex on K × R m .
By using Assumption 3 and the fact that

E(k) = 1 2 n i=1 f i (x(k), U (k)) 2
we can deduce that U j (k -s) → E(k) is convex for all 1 ≤ j ≤ m and 0 ≤ s ≤ r. So the backpropagation algorithm converges and

lim k→∞ ( ∂E(k) ∂U1(k-r) , . . . , ∂E(k) ∂Um(k-r) , . . . , ∂E(k) ∂U1(k) , . . . , ∂E(k) ∂Um(k) ) = 0.
Finally, Assumption 1 implies that lim k→∞ x(k) = 0. This can explain why the backpropagation algorithm is applied on the extended system (2) and not on the system (1). Thus, we have the following result:

Theorem 1. Under Assumptions 1-3, the solutions of the closed-loop system (1) with the DL feedback control given by the backpropagation algorithm ( 3)-(5) on the system (2) converge to the origin.

The DL feedback control approach is illustrated in Figure 2. are not available due to the fact that f is unknown, then we can use numerical differentiation methods for obtaining an approximate value of these derivatives [START_REF] Russell | Artificial Intelligence: A Modern Approach[END_REF]. Moreover, if r = 0, then Assumption 2 is replaced by the one that u j → f i (x, u) 2 is convex for all 1 ≤ i ≤ n and 1 ≤ j ≤ m.

If the system (1) is affine, i.e.,

x(k + 1) = φ(x(k)) + m j=1 g j (x(k))u j (k) (6)
with φ : R n → R n being continuous such that g j (x) = (g 1j (x), . . . , g nj (x)), then we have

∂E(k) ∂U j (k) = n i=1 x i (k + 1)g ij (x(k)) (7)
for all 1 ≤ j ≤ m. Note that affine systems are very common in applications [START_REF] Hasegawa | Control Problems of Discrete-Time Dynamical Systems[END_REF]. Finally, consider the very special case of linear systems of the form:

x(k + 1) = Ax(k) + Bu with A ∈ R n×n and B = [b ij ] ∈ R n×m , then we have ∂E(k) ∂U j (k) = n i=1 b ij x i (k + 1)
for all 1 ≤ j ≤ m. These results ensure that the backpropagation algorithms converge for the system (1) and the fundamental approximation theorem ensures that a DNN can approximate any continuous function on a compact set K. Now we can present how to implement the DL control algorithm.

Algorithm 1 explains how the backpropagation algorithm presented in Section III is used for building a DNN called DL feedback control. It stops when the draw of new random states in K does not make evolve the weights of the DNN. The learning process on the system (2) allows to find a DL feedback control for the system (1). If Algorithm 1 uses only the states {x(k) ∈ K : k ∈ N} of a solution of the system (1) without using random states, it will learn badly with too few states in K and therefore too little information. It is possible Algorithm 1: DL feedback control algorithm on K Input: x(k -r) ∈ K and > 0 Output: DNN U on K for the system (1) 1 Apply the backpropagation algorithm (3)-(5) to the system (2) with the input x(k -s) of the DNN where Remark 2. We use all the states (x(k), . . . , x(k -r)) in the input of the DNN for tuning its weights contrary to what is done in [START_REF] Li | Neural network control of unknown nonlinear systems[END_REF], [START_REF] Fourati | Stabilization of unknown nonlinear systems using neural networks[END_REF] where only both states (x(k), x(k -r)) are used under the restrictive assumption that the system (2) is uniquely invertible. Moreover an algorithm, namely Algorithm 1, is provided explaining how tuning the weights of the DNN which is not the case in [START_REF] Li | Neural network control of unknown nonlinear systems[END_REF], [START_REF] Fourati | Stabilization of unknown nonlinear systems using neural networks[END_REF].

0 ≤ s ≤ r. 2 if δw(k -s) 2 = i,j (k -s) 2 ≥ then 3 w ij (k -s) ← w ij (k -s + 1) 4 if E(k) < E(k -1) then 5 x(k -s) ← x(k -s + 1) for all 0 ≤ s ≤ r 6 if x(k -s + 1) ∈ K for all 0 ≤ s ≤ r then

IV. FORMATION CONTROL OF MULTI-AGENT SYSTEMS

In this part, we will present the formation control problem of multi-agent systems and show how to use Algorithm 1 to solve it. The formation control problem we consider here is based on the consensus protocol, which means each agent will have a copy of its own understanding consensus state x ci . However, the real state of each agent will differ from its understanding consensus state that the formation topology can be maintained. Here we consider the real state of the agent as x i (t) and each agent has its own understanding consensus state x ci . The relationship between x ci and x i can be written as:

x i (t) = x ci (t) + p i (8)
where p i is the vector denoting the real state of each agent x i (t) deviating its understanding consensus state x ci .

Moreover, the dynamic of the understanding consensus state can be written as follows:

x ci (t + 1) = x ci (t) + u i (t) (9)
for each agent. The goal of the consensus protocol is to reach

lim t→∞ x c1 (t) = • • • = lim t→∞ x cN (t) (10)
for all N agents. Next, we need to define the error dynamic for each agent i as:

e i (t) = x ci (t) -1/N N j=1 x cj (t) = 1/N N j=1 (x ci (t) -x cj (t)) (11)
Then, we have

e i (t + 1) = x ci (t + 1) -1/N N j=1 x cj (t + 1) = 1/N N j=1 (x ci (t) -x cj (t)) = 1/N N j=1 x ci (t) -x cj (t) + 1/N N j=1 u i (t) -u j (t) (12)
which means:

e i (t + 1) = e i (t) + 1/N N j=1 u i (t) -u j (t) (13)
Equation ( 13) can be written as:

e(t + 1) = e(t) + 1/N (I N -1 N 1 T N )u(t) (14)
where e(t) = (e 1 (t), . . . , e N (t)) T , u(t) = (u 1 (t), . . . , u N (t)) T , 1 N is a column vector witch has N elements with each element equals 1. The system error ( 14) can be seen as the specific form of the dynamical system (2) with the function f = e(t) + 1/N (I N -1 N 1 T N )u(t). Thus, our goal has become to design a DNN controller based on Algorithm 1 to force the error state e(t) to converge to 0.

V. EXPERIMENT

In this part, we will implement an experiment to show the effectiveness of the Algorithm 1 on our multirobot platform. We first use the Gazebo simulation environment to train the DNN controller to obtain the network model. After that we could use the trained network in the real robot to show the effectiveness of Algorithm 1. The multirobot system used here is 4 Pioneer 3dx robots. Each robot is with a Nvidia TX1 so that they can communicate with each other by the Robot Operating System (ROS). The DNN we used here has 6 hidden layers and each layer is composed of 100 neurons. We choose the ReLU function as the activation function. After training the network, the robot system can achieve the formation control by walking along the path generated by the algorithm.

Let (r xi , r yi ), θ i denote the Cartesian position and orientation of the ith robot, respectively. Thus, we can denote the state of the robot x i = [r xi , r yi , θ i ] T . The linear and angular speed of the robot can be denoted by (v i , w i ). Thus, we have the kinematic equation of the ith robot as follows:

ṙxi = v i cos(θ i ), ṙyi = v i sin(θ i ), θi = w i (15)
Then, by linearizing (15) around a fixed point off the center of the wheel axis (x i , ȳi ) of the robot, where xi = r xi + d i cos(θ i ), ȳi = r yi + d i sin(θ i ), and d = 0.15m, we have the equation

v i w i = cos(θ i ) sin(θ i ) -(1/d) sin(θ i ) (1/d) cos(θ i ) a xi a yi .
This equation is a simple kinematic one that can transfer the action of each robot to the velocity command. Next, we can define the understanding consensus state for each robot as x ci = [x r ci , y r ci , θ r ci ] T . The [x r ci , y r ci ] and θ r ci denote the reference position and orientation of the understanding formation center of the robots team for each robot i, respectively. Since the formation center of the robots team is changing dynamically, each robot would maintain a local variable x ci , which is the sensing value of the state s r by each robot i.

Moreover, there is a relationship (8) that we should clarify between the real state x i and the understanding consensus state x ci :

xi ȳi = x r ci y r ci + cos(θ r ci ) -sin(θ r ci ) sin(θ r ci ) cos(θ r ci ) xif ỹif , (16)
where [x if , ỹif ] T represents the desired deviation vector of the ith robot relative to the geometric center of the formation. The initial positions for the robots are x1 = 0.935m, y1 = 0.935m, x2 = 2.205m, y2 = 0m, x3 = 1.46m, y3 = 1.46m, and x4 = -0.81m, y4 = 2.32m. Also, we add a virtual leader for the robots team and its initial position is x0 = 0m and y0 = 0m. Its velocity is v 0x = 0.01m/s and v 0y = 0.01m/s. Also, x1f = 1.5m, ỹ1f = 1.5m, x2f = -1.5m, ỹ2f = 1.5m, x3f = 1.5m, ỹ3f = -1.5m, and x4f = -1.5m, ỹ4f = -1.5m. Figures 3 and4 show the trajectory of each robot moving in x and y direction, respectively. Moreover, Figures 5 and6 show the velocity of each robot in x and y direction, respectively. Next, we will perform the experiment on the real robot by using the trained DNN before. Figure 7 shows the robots are in the initial position, Figure 8 shows that the robots are running, and Figure 9 shows that they have achieved the formation finally.

VI. CONCLUSION

In this paper, we consider the formation control problem by using a deep learning based feedback control method. The backpropagation algorithm is used to train the DNN to obtain the controller. In order to implement this method, we transfer the formation control problem into a general problem (2) that can be solved by Algorithm 1. Finally, we verified the effectiveness of the algorithm on our multirobot platform.

Fig. 1 .

 1 Fig. 1. Classical feedback control.

Fig. 2 .

 2 Fig. 2. DL feedback control.

  convergence speed of the DL feedback control by replacing Step 4 of Algorithm 1 with E(k) < λE(k -1) where 0 < λ < 1.

Fig. 3 .

 3 Fig. 3. The x direction trajectory of each robot.

Fig. 4 .

 4 Fig. 4. The y direction trajectory of each robot.

Fig. 5 .

 5 Fig. 5. The velocity of x direction trajectory of each robot.

Fig. 6 .

 6 Fig. 6. The velocity of y direction trajectory of each robot.