A Nicolas
email: nicolas@geologie.ens.fr

G Blöcher

C Kluge

Z Li

H Hofmann

L Pei

H Milsch

J Fortin

Y Guéguen

Pore pressure pulse migration in microcracked 1 andesite recorded with fibre optic sensors 2

Pore pressure has a major influence on the effective stress and thus on the mechanical behaviour of rocks. In this study, we focus on the hydromechanical behaviour of a low porosity andesitic rock heat-treated to 930 o C to induce thermal cracks and increase the permeability of the samples. First, we show that permeability decreases from 8 × 10 -16 m 2 to 1.5 × 10 -17 m 2 with a confining pressure (P c ) increase from 2 MPa to 40 MPa (pore pressure being approximately 0.2 MPa). Then, we used fibre optic pressure sensors to monitor pore pressure diffusion at three points along the sample during the propagation of a pore pressure pulse under hydrostatic (P c =40 MPa) and triaxial stresses (P c =40 MPa, differential stress of 356 MPa). When the pore pressure pulse was applied, the fibre optic sensors showed a sudden pore pressure increase one after the other as a function of their location along the sample. Pore pressure increase downstream was very smooth under hydrostatic stress and almost zero after the duration of the experiment (50 minutes) under triaxial stresses. This lack of downstream pore pressure increase under triaxial stresses is due to the fact that a differential stress of 356 MPa decreased permeability from approximately 10 -17 m 2 to approximately 10 -19 m 2 . Finally, the pore pressure diffusion process was modelled considering a uniform spatial distribution of permeability in the andesite sample and the dead volume attached at the downstream side.

Introduction

Rock permeability and its evolution with changes of external stress and pore pressure is of interest for industrial purposes such as petroleum and geothermal reservoir engineering, radioactive waste management, CO 2 sequestration and for the understanding of natural hazards such as earthquakes. Human activities can lead to increases in effective stress because of pore pressure decrease, for example during hydrocarbon or geothermal fluid extraction [e.g. [START_REF] Segall | A note on induced stress changes in hydrocarbon and geothermal reservoirs[END_REF]]. Increases of effective pressure (external stress minus pore pressure) can lead to a decrease of permeability [e.g. [START_REF] David | Laboratory measurement of compaction-induced permeability change in porous rocks: Implications for the genera-D R A F T January[END_REF][START_REF] Miller | Properties of large ruptures and the dynamical influence of fluids on earthquakes and faulting[END_REF][START_REF] Bemer | From injectivity to integrity studies of CO 2 geological storage chemical alteration effects on carbonates petrophysical and geomechanical properties[END_REF] causing a reduction of the production rate in geo-fluids exploitation. In geothermal fields situated in andesitic areas, permeability distribution can vary over a wide range of values and fault zones usually control fluids pathways [e.g. Brehme et al., 2016a, b]. As geothermal energy is currently developed [e.g. [START_REF] Tester | The future of geothermal energy: Impact of enhanced geothermal systems (EGS) on the United States in the 21st century[END_REF][START_REF] Zimmermann | Hydraulic stimulation of a deep sandstone reservoir to develop an enhanced geothermal system: Laboratory and field experiments[END_REF], it is necessary to understand the links between fluid pressure, mechanical behaviour and possible seismicity at all scales [e.g. Shapiro and Dinske, 2009a, b;[START_REF] Zang | Analysis of induced seismicity in geothermal reservoirs-an overview[END_REF][START_REF] Johann | Scaling of seismicity induced by nonlinear fluid-rock interaction after an injection stop[END_REF][START_REF] Katayama | Fluid-induced fracturing of initially damaged granite triggered by pore pressure buildup[END_REF][START_REF] Chanard | Sensitivity of acoustic emission triggering to small pore pressure cycling perturbations during brittle creep[END_REF]. In enhanced geothermal systems, stimulation significantly increases pore pressure [START_REF] Häring | Characterisation of the Basel 1 enhanced geothermal system[END_REF], leading to a reduction of effective stress. Reduced effective stress generates new fluid path ways, thus increasing production efficiency [START_REF] Zimmermann | Pressure-dependent production efficiency of an enhanced geothermal system (EGS): stimulation results and implications for hydraulic fracture treatments[END_REF] but it may also induce slip on pre-stressed discontinuities [e.g. [START_REF] Majer | Induced seismicity associated with enhanced geothermal systems[END_REF][START_REF] Deichmann | Earthquakes induced by the stimulation of an enhanced geothermal system below Basel (Switzerland)[END_REF][START_REF] Terakawa | High fluid pressure and triggered earthquakes in the enhanced geothermal system in Basel, Switzerland[END_REF].

In the laboratory, it has been shown that even in the elastic domain, the effective stress law for the permeability is not straightforward [e.g. [START_REF] Ghabezloo | Effective stress law for the permeability of a limestone[END_REF][START_REF] Braun | Theoretical analysis of pore pressure diffusion in some basic rock mechanics experiments[END_REF] but is of major importance for prediction of permeability evolution in stressed geo-materials [e.g. [START_REF] Shapiro | Estimating the crust permeability from fluid-injection-induced seismic emission at the KTB site[END_REF][START_REF] Shapiro | Triggering of seismicity by pore-pressure perturbations: Permeability-related signatures of the phenomenon[END_REF]. At stresses beyond the elastic domain of X -4 NICOLAS ET AL.: FIBRE OPTIC LOCAL PORE PRESSURE the rock, crack propagation or nucleation, and plasticity add complexities. Depending on the failure mode of the rock considered, permeability can either increase or decrease with increasing damage [e.g. [START_REF] Zhu | The transition from brittle faulting to cataclastic flow: Permeability evolution[END_REF] because rock deformation can occur with compaction or dilatancy of the porous medium. Brittle deformation can lead to dramatic changes of permeability because of changes of crack density and connectivity [e.g. [START_REF] Guéguen | Models and time constants for permeability evolution[END_REF][START_REF] Darot | Acoustic wave velocity and permeability evolution during pressure cycles on a thermally cracked granite[END_REF]Guéguen et al., 2011;[START_REF] Nicolas | Radon emanation from brittle fracturing in granites under upper crustal conditions[END_REF]. The sustainability of fracture permeability under varying pressure is questionable [e.g. Zimmermann et al., 2010;[START_REF] Kluge | Sustainability of fractured rock permeability under varying pressure[END_REF]. When compaction is taking place, permeability can decrease with increasing deformation and damage [e.g. [START_REF] Fortin | Elastic wave velocities and permeability evolution during compaction of Bleurswiller sandstone[END_REF][START_REF] Farquharson | Inelastic compaction and permeability evolution in volcanic rock[END_REF].

Permeability evolution as a function of stress and strain in volcanic rocks has already been studied and can decrease [e.g. [START_REF] Farquharson | Inelastic compaction and permeability evolution in volcanic rock[END_REF] or increase [e.g. Farquharson et al., 2016b] with increasing axial strain. Moreover, in active volcanic hydrothermal systems, rocks can have very different petrophysical properties. For example, [START_REF] Heap | A multidisci-D R A F T January[END_REF] showed that the porosity of the materials in an active volcanic hydrothermal system (in this study the Whakaari stratovolcano) varies from 0.01 to 0.7 and permeability varies by eight orders of magnitude (from 10 -19 m 2 to 10 -11 m 2 ). The authors assumed that the wide range in physical and hydraulic properties is the result of the numerous lithologies and their varied micro-structures and alteration intensities. Finally, some post-earthquake ground movement might be correlated to pore-pressure transients [e.g. [START_REF] Jonsson | Post-earthquake ground movements correlated to pore-pressure transients[END_REF], implying that pore pressure pulse diffusion under varying permeability conditions is of interest.

Fibre optic pressure sensors are used in boreholes to record pressure with a high sensitivity [e.g. [START_REF] Reinsch | Temperature-dependent characterization of optical fibres for distributed temperature sensing in hot geothermal wells[END_REF]Reinsch et al., 2013a[START_REF] Reinsch | On the mechanical coupling of a fiber optic cable used for distributed acoustic/vibration sensing applicationsa theoretical consideration[END_REF]. Fibre optic Fabry-Perot interferometer diaphragm and the method of measurement is known since the early 20th D R A F T January 20, 2020, 9:18pm D R A F T X -5

century, but to the authors' knowledge, were not widely used in laboratory experiments [e.g. Reinsch et al., 2012a;[START_REF] Blöcher | Direct and indirect laboratory measurements of poroelastic properties of two consolidated sandstones[END_REF].

In the frame of the GEOTREF project, aiming at improving our understanding of the behaviour of fractured geothermal reservoirs, this study focuses on the hydro-mechanical behaviour of a low porosity andesite rock. Samples were heat-treated to 930 o C to induce thermal cracks and increase their permeability. The questions we address are: How does a pulse of pore pressure propagate in a low porosity and low permeability rock? What are the induced spatial distribution and evolution of the permeability under hydrostatic and triaxial conditions? Can fibre optic Fabry-Perot interferometer diaphragms be used in laboratory rock deformation experiments to record local pore pressure evolution and nonlinear diffusion of pore pressure? We report results of local pore pressure and permeability measurements during a pore pressure pulse propagation under hydrostatic (P c =40 MPa) and triaxial stress states (P c =40 MPa, differential stress of 356 MPa).

Sample characterization and experimental apparatus

Rock material and sample preparation

A cylindrical sample (length=100mm, diameter=50mm) was cored from a block of andesite recovered from an outcrop in Anse a la Barque (Guadeloupe, French West Indies)

and the end-faces were ground flat to ensure a good parallelism. This andesite is formed of a groundmass (composed of small crystals) with large phenocrysts (from hundreds of microns to a few millimetres for the largest plagioclases) commonly micro-cracked (Figure 1 A and B). No obvious bedding related to volcanic features is identified in the block. The sample has a connected porosity of less than 1% as measured using a pulse method. Initial P-wave velocity in the sample was measured at room pressure and temperature and was around 5000 m/s. Permeability is in the order of 10 -21 m 2 as measured using a constant flow procedure [START_REF] Li | Physical and mechanical properties of thermally cracked andesite under pressure[END_REF]. To induce thermal cracks and increase its permeability before the experiment, the sample was heat-treated with a programmable Meker MHT-3 furnace [START_REF] Nicolas | Radon emanation from brittle fracturing in granites under upper crustal conditions[END_REF]. The sample was heated at a rate of 50 o C/h, up to 930 o C, kept 2 h at this temperature, and then cooled at 50 o C/h, to avoid any quenching or thermal shock effect. Details on this procedure applied on andesite and associated discussion can be found in [START_REF] Li | Physical and mechanical properties of thermally cracked andesite under pressure[END_REF]. The sample shows a higher fracture density after heat treatment at 930 o C compared with intact samples (Figure 1). Heat treatment induces a decrease of P-wave velocities to approximately 2900 m/s and an increase of permeability to 10 -15 m 2 measured using a constant flow procedure [START_REF] Li | Physical and mechanical properties of thermally cracked andesite under pressure[END_REF], as also observed by [START_REF] Darot | Permeability of thermally cracked granite[END_REF] among others. The intense micro-fracturing is due to the different thermal expansion coefficients of the minerals [START_REF] Fredrich | Micromechanics of thermally induced cracking in three crustal rocks[END_REF][START_REF] Browning | Cooling-dominated cracking in thermally stressed volcanic rocks[END_REF].

Experimental apparatus and fibre optic sensor recording system

The sample was deformed in a conventional triaxial testing device (MTS 815) installed at the GFZ Postdam (Germany). Details on the triaxial apparatus can be found in [START_REF] Blöcher | Investigation of the undrained poroelastic response of sandstones to confining pressure via laboratory experiment, numerical simulation and analytical calculation[END_REF]. A shrinking tube made of fluorinated ethylene propylene (FEP)

was used during measurements of permeability evolution under hydrostatic loading and neoprene jackets were used during pore pressure pulse propagation experiments to separate the sample from the oil confining medium. Axial and radial strains were recorded with two axial extensometers and one radial chain extensometer directly attached to the jacket. Axial displacement was also measured with displacement transducers (DCDT) mounted between the moving piston and a fixed platen. Volumetric strain is calculated as ε v = ε ax + 2ε r where ε ax and ε r are the axial and radial strains, respectively. From now on, compressive stresses and compactive strains are defined as positive.

Besides the mechanical system, an independent hydraulic system dedicated to pore pressure was connected to the sample. Pore pressure and fluid flow can be controlled at the top

D R A F T January 20, 2020, 9:18pm D R A F T NICOLAS ET AL.: FIBRE OPTIC LOCAL PORE PRESSURE X -7
and the bottom of the sample with two servo-controlled micro-volumetric pumps (Quizix 6000-Series) coupled by pairs with an accuracy of about 10 -2 MPa and 10 -4 ml/min, respectively [e.g. [START_REF] Blöcher | Direct and indirect laboratory measurements of poroelastic properties of two consolidated sandstones[END_REF]. The total capillary void space between the pumps and the sample is approximately 100 ml. The triaxial cell is equipped with a heating system. The temperature inside the vessel is recorded by two thermo-couples plunged in the confining oil close to the sample [e.g. [START_REF] Pei | Thermomechanical properties of Upper Jurassic (Malm) carbonate rock under drained conditions[END_REF].

The fibre optic sensors used in this study consist of miniature all-silica extrinsic Fabry-Perot cavity interferometer (EFPI) pressure sensors. Monochromatic light is sent in the fibre and propagates to its head. Incident light is reflected two times: (1) it is first partly reflected at the entrance into the EFPI cavity (glass/air interface) and (2) the transmitted light becomes reflected at the termination of the cavity (air/glass interface). Light reflected from the end of the cavity is partly transmitted back into the fibre and interferes with light reflected from the first reflection at the entrance into the Fabry-Perot cavity.

The phase shift between the two reflected signals depends on the cavity length, which is a function of pressure and temperature, the latter being maintained constant at 30 o C during all the experiments. Before being used, the fibre optic sensors have to be calibrated (see next section), that is to say that the relation between pressure and cavity length (and thus the phase shift evolution with pressure) has to be known through a calibration. A thorough theoretical discussion of the principle of this sensor can be found in [START_REF] Bremer | Fibre optic pressure and temperature sensor for geothermal wells[END_REF] and details on the experimental set-up are given in Reinsch et al. [2012a, b].

Feedthroughs constructed in-home were specifically dedicated to let fibre optics enter the pressure vessel [Reinsch et al., 2013b]. To measure the pore pressure in the rock specimen submitted to confining pressure, three fibre optic sensors were placed into three holes that were previously drilled through the neoprene jacket in the rock at 1/4, 2/4 and 3/4 of the height of the sample (Figure 2). The combination of the fibre optic sensors with the D R A F T January 20, 2020, 9:18pm D R A F T pressure sensors of the Quizix system allow five independent pressure measurements at five points along the sample (0, 25, 50, 75 and 100 mm).

Calibration of the fibre optic sensors

Before using the fibre optic sensors, these sensors had to be calibrated. For calibration, the three fibre optic sensors were placed within the oil in the pressure chamber. Temperature of the oil was controlled and maintained constant at 30 o C. The pressure was increased from 2 to 40 MPa at a rate of 0.25 MPa/min in the vessel while data from the sensors was recorded with a temporal resolution of 15 s, that is to say one measurement approximately each 0.06 MPa. The sensor responses to an increase of pressure were tested and the sensors response is shown to be linear as a function of pressure (Figure 3) for two of the sensors (Fibres 1 and 3). Fibre 2 shows a hysteresis that could be explained by a small leak in the Fabry-Perot cavity. The response of Fibre 2 seems to be correct and linear after an initial non-linear response (Figure 3). From now on, this sensor will be used only in the linear response domain and all pressure calculations are derived by the slope of the calibration.

With the calibration information, it is possible to correlate phase shift changes to pressure changes. Calibration information was then used to infer the pressure evolution as a function of the sensor responses during the experiments. Based on the hysteresis of fibres 1 and 3, the precision can be estimated to be around 0.5 MPa. It is more complicated to estimate the precision of fibre 2, and the response of this fibre will always be shown as a dashed line.

Experimental procedures

In this section, we outline the experimental procedures used to monitor (1) the evolution of porosity and permeability during hydrostatic loading (section 3.1), (2) the pore pressure

D R A F T January 20, 2020, 9:18pm D R A F T NICOLAS ET AL.: FIBRE OPTIC LOCAL PORE PRESSURE X -9
evolution while a pore pressure pulse is sent under hydrostatic pressure (section 3.2), and

(3) the pore pressure evolution while a pore pressure pulse is sent under triaxial stresses (section 3.3). One single sample was used for the three different phases of the experiment.

Temperature of the oil was controlled and maintained constant at 30 o C during all the phases of the experiments.

3.1. Measuring the evolution of porosity and permeability during hydrostatic loading 3.1.1. Experimental procedure for hydrostatic loading

The dry sample was first placed in the chamber and vacuumed to approximately 10 mbar. It was then loaded hydrostatically up to a pressure of 2 MPa and saturated with deionized water. Saturation was controlled with the micro-volumetric pumps. Full saturation was assumed to be reached when the total injected volume of water remained constant. The sample was loaded hydrostatically up to a confining pressure of 40 MPa at a controlled pressure increase rate of approximately 0.25 MPa/min, which leads to a maximum volumetric strain rate of approximately 4.4 × 10 -7 s -1 if we consider a bulk modulus of K = 9.5 GPa (as measured, see section 4.1). The characteristic time t to reach fluid pressure equilibrium for a diffusion process over a distance l can be approximated by [START_REF] Carslaw | Conduction of heat in solids[END_REF][START_REF] Ge | Hydrodynamic response to strike-and dip-slip faulting in a half-space[END_REF][START_REF] Duda | The weakening effect of water on the brittle failure strength of sandstone[END_REF][START_REF] Nicolas | Brittle and semibrittle behaviours of a carbonate rock: influence of water and temperature[END_REF]:

t ∼ l 2 D , ( 1
)
where D is the hydraulic diffusivity. The hydraulic diffusivity D can be calculated as [START_REF] Kümpel | Poroelasticity: parameters reviewed[END_REF]:

D ∼ kBK d ηα , ( 2
)
where k is the permeability, B is Skempton's coefficient, K d is the drained bulk modulus, η is the fluid viscosity, and α is Biot's coefficient. For heat-treated andesite saturated

D R A F T January 20, 2020, 9:18pm D R A F T X -10 NICOLAS ET AL.: FIBRE OPTIC LOCAL PORE PRESSURE
with water, k = 10 -17 m 2 (lowest value of permeability measured during hydrostatic loading, see section 4.1), η = 1.002 × 10 -3 Pa.s, and K d = 9.5 GPa. Assuming that Skempton's coefficient and Biot's coefficient are of the order of unity, hydraulic diffusivity is D ∼ 9.5 10 -7 m 2 s -1 . The length of the sample is l = 10 cm, and thus the diffusion time is in the order of t ∼ 100 s, which is a shorter time than hydrostatic loading time (t ∼ 10000 s) and ensures a fully drained behaviour.

Once the confining pressure of 40 MPa was reached, the sample was hydrostatically unloaded using the same procedure as that used for loading. The sample is then removed from the pressure vessel and three holes are drilled to place the three fibre optic sensors (see section 2.2).

Permeability measurements during hydrostatic loading

During hydrostatic loading, permeability is continuously recorded. To do so, the hydraulic short-cut (capillary connection) between the two ends of the sample were hydraulically disconnected. Pore pressure was kept at a constant pressure of 0.2 MPa on the downstream side and a constant flow was applied on the upstream side of the sample.

First a constant flow of 0.4 ml/min was used, and was then decreased to 0.05 ml/min, 0.025 ml/min and finally 0.01 ml/min, as permeability decreased and to prevent upstream pressure to overcome confining pressure. A constant flow of 0.01 ml/min was used during complete unloading. Pore pressure increased on the side at which the flow was applied and permeability was calculated using Darcy's law [START_REF] Darcy | Les fontaines publiques de la ville de Dijon: exposition et application[END_REF]:

k = QηL A∆P , ( 3
)
where k is the permeability, Q is the flow through the sample, η is the fluid viscosity, L is the length of the sample, A is the sample cross-sectional area and ∆P is the pore pressure differential.

3.1.3. Can compaction induce a flow comparable to the imposed one?

D R A F T January 20, 2020, 9:18pm D R A F T NICOLAS ET AL.: FIBRE OPTIC LOCAL PORE PRESSURE X -11
Using the experimental procedure described above, a remaining question is whether compaction could lead to a flow comparable to the one imposed by the permeability measurement [e.g. [START_REF] Faulkner | The effect of temperature, the nature of the pore fluid, and subyield differential stress on the permeability of phyllosilicate-rich fault gouge[END_REF][START_REF] Mitchell | Experimental measurements of permeability evolution during triaxial compression of initially intact crystalline rocks and implications for fluid flow in fault zones[END_REF][START_REF] Crawford | Strength, porosity, and permeability development during hydrostatic and shear loading of synthetic quartz-clay fault gouge[END_REF][START_REF] Fortin | Influence of thermal and mechanical cracks on permeability and elastic wave velocities in a basalt from Mt. Etna volcano subjected to elevated pressure[END_REF]. Volumetric strain rate is calculated as a function of imposed hydrostatic stress rate and bulk modulus as follows:

εv = Ṗ K , ( 4
)
where Ṗ is the controlled pressure increase rate. Compactive volume rate is then simply calculated as:

∆ V = V εv , ( 5
)
where V is the volume of the sample. Using V = 200ml and εv = 2.6 × 10 -5 min -1 , it yields an approximate compactive volume rate of ∆ V = 5 × 10 -3 ml/min. This approximate (upper bound) compactive volume rate is twice lower than the imposed flow for the permeability measurement in the sample. Assuming a Biot's coefficient close to unity will forward the bulk volume rate to the fluid volume rate. For higher loads, Biot's coefficient will decrease, leading to a reduction of the fluid volume rate induced by compaction.

To confirm that the compactive volume rate is lower than the imposed flow for the per- The pore pressure pulse propagation experiment was performed at a confining pressure of 40 MPa. Once the fibre optic sensors were plugged in the previously 25 mm deep drilled holes, the dry sample was first placed in the pressure vessel and loaded isostatically up to a pressure of 2 MPa. Vacuum was then applied in the total pore pressure system. Once a pore pressure of 0.01 bar was attained, the sample was saturated with deionized water.

Pore pressure was then maintained constant at 0.2 MPa on both sides of the sample.

Saturation was controlled with the micro-volumetric pumps. Full saturation was assumed to be reached when the total injected volume of water remained constant.

Once full saturation was reached, pore pressure was maintained constant at a pressure of 0.2 MPa on both sides of the sample and the hydrostatic pressure was increased from 2

MPa to 40 MPa at a controlled pressure increase rate of approximately 0.25 MPa/min, ensuring a fully drained behaviour (see section 3.1.1). When the total injected volume of pore water remained constant, the hydraulic short-cut (capillary by-pass) between the two ends of the sample were hydraulically disconnected and pore pressure was increased from 0.2 MPa to 20 MPa (in a few seconds) and kept constant at 20 MPa on the bottom side of the sample (from now on called upstream) while the other pump (from now on called downstream) was stopped and did not regulate the pore pressure. Downstream pump was used to record the evolution of the pore pressure on this side of the sample. Pore pressure evolution was recorded on both sides of the sample with the Quizix pumps and on three points along the sample with the fibre optic sensors, leading to five independent pressure measurements along the sample (0, 25, 50, 75 and 100 mm).

Once pore pressure reaches equilibrium in the whole sample, pore pressure was decreased and maintained constant at a pressure of 0.2 MPa on the downstream side of the sample while the pump situated on the upstream side was stopped. Finally, 35 minutes later,

D R A F T January 20, 2020, 9:18pm D R A F T X -13
pore pressure was decreased and maintained constant at a pressure of 0.2 MPa on both sides of the sample.

Pore pressure pulse under triaxial stresses

First, the saturated sample was maintained under a hydrostatic pressure of 40 MPa following the procedure described in section 3.2. Then, differential stress was applied by axially loading the sample at a controlled stress rate of ∼ 10 -1 MPas -1 until a differential stress of 356 MPa was attained. This controlled stress rate ensures a fully drained behaviour during loading. Under the target differential stress, the sample remains in the elastic domain as shown by the results released in [START_REF] Li | Physical and mechanical properties of thermally cracked andesite under pressure[END_REF]. Once the target differential stress was reached, the stress-state was maintained constant. The two ends of the sample were hydraulically disconnected and a pressure pulse was then induced. Pore pressure was increased from 0.2 MPa to 35 MPa (in a few seconds) and kept constant at 35 MPa on the upstream side of the sample while the downstream pump was stopped.

Again, the pore pressure evolution was recorded on both sides of the sample with the Quizix pumps and on three points along the sample with the fibre optic sensors, leading to five independent pressure measurements along the sample (0, 25, 50, 75 and 100 mm).

Results

Hydrostatic loading: compaction and permeability evolution

Results for the hydrostatic loading and unloading experiment with up to 40 MPa confining pressure are presented in Figure 5 (a). During loading, the hydrostatic response was non-linear up to a pressure of ∼ 28 MPa, beyond which the stress-strain curve became almost linear with a slope corresponding to a static bulk modulus of approximately K = 9.5 GPa (Figure 5 a). The mechanical response of heat-treated andesite subjected to isostatic loading is typical of a rock with micro-cracks [START_REF] Walsh | The effect of cracks on the compressibility of rock[END_REF]. The non-linearity [START_REF] Walsh | The effect of cracks on the compressibility of rock[END_REF][START_REF] Batzle | Microcrack closure in rocks under stress: Direct observation[END_REF][START_REF] Baud | Dilatancy, compaction, and failure mode in Solnhofen limestone[END_REF][START_REF] Vajdova | Compaction, dilatancy, and failure in porous carbonate rocks[END_REF][START_REF] Nicolas | Brittle and semibrittle behaviours of a carbonate rock: influence of water and temperature[END_REF], such as the ones shown in Figure 1 (B, C and D). All cracks are considered to be closed at an isostatic stress above 28 MPa. Comparing the volumetric strain measured and the perfectly elastic trend, micro-crack porosity can be estimated to be at least mp ∼ 0.17% [START_REF] Walsh | The effect of cracks on the compressibility of rock[END_REF]. During unloading, the hydrostatic response was non-linear and volumetric strain did not reach 0 when pressure reached 2 MPa.

D R A F T January 20, 2020, 9:18pm D R A F T
As pressure was increased from 2 to 40 MPa, permeability decreased from 8 × 10 -16 m 2 to 1.5 × 10 -17 m 2 (Figure 5 b). Together with the static measurements, this decrease highlights the closure of pre-existing cracks [e.g. [START_REF] Brace | Permeability of granite under high pressure[END_REF][START_REF] Brace | Permeability from resistivity and pore shape[END_REF][START_REF] Brace | A note on permeability changes in geologic material due to stress[END_REF] between 0 and 40 MPa and more especially between 2 MPa and 28 MPa when the decrease is more dramatic. During unloading, permeability increased from its minimum value to 2 × 10 -16 m 2 (Figure 5 b) but did not recover its initial value. This hysteresis, together with the mechanical data, implies that some micro-cracks remain closed after the loadingunloading procedure.

The evolution of permeability as a function of volumetric strain is strictly identical during loading and unloading (Figure 5 c), showing that the cracks that remaining closed after unloading are responsible for the permanent decrease of permeability after unloading.

Pore pressure propagation under hydrostatic stress

Local pore pressure measurements with fibre optic sensors as a function of time are presented in Figure 6 (a). In addition, the evolution of upstream and downstream pore pressures in the pumps are plotted in blue and black, respectively. During the whole pressure pulse propagation, hydrostatic pressure was 40 MPa. Then, pore pressure was decreased suddenly at the former downstream point (Figure 6 (a)), leading to a negative pore pressure pulse symmetrical to the positive pulse sent at the beginning of the experiment. Pore pressure decreased and similar to the first part of the experiment, a pore pressure gradient was highlighted. After 35 minutes, one can observe a small positive pore pressure pulse during global pore pressure decrease, which is due to the opening of a valve outside the sample.

After 65 minutes, pore pressure was maintained at 0.2 MPa on both sides of the sample.

The fibre optic sensors showed a progressive decrease of pore pressure in the sample. The middle fibre optic sensor did not decrease to zero, which could be due to its hysteresis shown during the calibration (section 2.3). that axial and volumetric strains were extremely low (tending to 0) implying that there was no creep.

Pore pressure propagation under triaxial stresses

When pore pressure was increased upstream, the fibre optic sensors reacted one after the other as a function of their location along the sample but downstream pressure remained at 0 MPa (Figure 8). The fibre optic sensors first showed a sudden increase to intermediate values ranging from 7.5 MPa (the nearest from downstream) to approximately 27.5 MPa (nearest to upstream). Pore pressure increase downstream was almost zero because pore pressure at the upstream side was decreased before it started increasing downstream (Fig-

ure 8). Pore pressures then remained constant during more than 50 minutes and between the two extreme (downstream and upstream) pore pressures. The gradient of pressure was approximately constant within the sample although effective pressure increased from the upstream side to the downstream side. After 55 minutes, pore pressure was maintained at 0.2 MPa at the former upstream point and pore pressure started to decrease all over the sample (Figure 8). The behaviour was similar as that observed under hydrostatic pressure, with the sudden decrease followed by a gradual decrease.

Discussion

Permeability evolution under hydrostatic stress

The mechanical response of the heat-treated andesite subjected to isostatic loading (Figure 5 a) is typical for a rock with micro-cracks and pores [e.g. [START_REF] Walsh | The effect of cracks on the compressibility of rock[END_REF], in agreement with micro-structural observations of the intact heat-treated rock (Figure 1). The At crack closure, permeability was about 1.5 × 10 -17 m 2 , whereas the thermally untreated sample had a permeability of approximately 10 -21 m 2 . At crack closure, the difference to the thermally untreated sample remained around four orders of magnitude, as also observed by [START_REF] Nara | Influence of macro-fractures and micro-fractures on permeability and elastic wave velocities in basalt at elevated pressure[END_REF]. This could be explained by the fact that cracks may be "mechanically closed", that is to say that asperities on both sides of the crack are in contact but there remains some open space between the asperities as observed by [START_REF] Kranz | The permeability of whole and jointed Barre granite[END_REF]; [START_REF] Bernabe | The effective pressure law for permeability in Chelmsford granite and Barre granite[END_REF] and modelled by [START_REF] Walsh | Effect of pore pressure and confining pressure on fracture permeability[END_REF]; [START_REF] Zimmerman | The effect of contact area on the permeability of fractures[END_REF]. [START_REF] Sevostianov | On elastic compliances of irregularly shaped cracks[END_REF] showed that asperities in contact (also called islands)

dominate the mechanical response.

To interpret theoretically and quantitatively the decrease of volumetric strain and permeability, we consider a porous rock as made of a mixture of solid matrix, spherical pores and penny-shaped cracks. Following [START_REF] Walsh | The effect of cracks on the compressibility of rock[END_REF], the crack closure pressure P cl for an isotropic stress state can be related to the crack aspect ratio defined as ξ = w/2c, where w and 2c are the crack aperture and the crack length, respectively :

ξ = 4P cl (1 -ν 2 0 ) πE 0 , (6)
where E 0 and ν 0 are the bulk moduli and Poisson's ratio of the crack-free matrix, respectively. Using the elastic properties of the crack-free material (E 0 = 50 GPa and ν 0 = 0.25, which are standard values) and a crack closure pressure of P cl = 28 MPa, equation 6 gives an aspect ratio equal to ξ = 7 × 10 -4 , in good agreement with values reported in the literature [e.g. [START_REF] Fontaine | Influence of melt viscosity of basaltic and andesitic composition on seismic attenuation in partially molten gabbronorite[END_REF][START_REF] Adelinet | Frequency and fluid effects on elastic properties of basalt: Experimental investigations[END_REF]Guéguen et al., 2011]. [START_REF] Adelinet | Deriving microstructure and fluid state within the Icelandic crust from the inversion of tomography data[END_REF] showed that at larger scale, aspect ratios are higher, highlighting that the aspect ratios could be scale-dependent in volcanic rocks. If cracks are assumed to be characterized by a penny-shape geometry, then crack porosity is given by [START_REF] Guéguen | Effective elastic properties of cracked rocks -An overview[END_REF]:

φ crack = 2πρ c ξ , (7)
D R A F T January 20, 2020, 9:18pm D R A F T

X -18 NICOLAS ET AL.: FIBRE OPTIC LOCAL PORE PRESSURE

where ξ is the average crack aspect ratio defined as ξ = w/2c and ρ c is the crack density defined as:

ρ c = 1 V n i=1 c 3 i , (8)
where c i is the radius of the ith crack and N is the total number of cracks embedded in the representative elementary volume (REV) V . Taking ξ = 7 × 10 -4 and an initial crack porosity mp 0.17% (determined from hydrostatic loading mechanical data shown in Figure 5 a), one finds an initial crack density equal to ρ c = 0.39.

Dienes et al.

[1982] calculated the permeability of a medium containing an isotropic distribution of cracks of radius c and found that:

k = 32π 15 ξ 3 c 2 ρ c θf, (9)
where θ and f are accounting for the hydrodynamics of flow through a system of cracks with varying thickness and accounting for the fraction of cracks that belong to an infinite network and will be determined from percolation theory, respectively. The percolation threshold calculated by [START_REF] Guéguen | Transport properties of rocks from statistics and percolation[END_REF] for the case of a Bethe lattice where each crack has 4 neighbours is a crack density of approximately 0.14, much lower than the crack density of 0.39 calculated for our heat-treated andesite. Above the percolation threshold, the function f is equal to 1, which is the value used here. The factor θ accounting for the hydrodynamics of flow through a system of cracks with varying thickness is taken equal to 1 here, in agreement with [START_REF] Guéguen | Transport properties of rocks from statistics and percolation[END_REF]. Using equation 9

with k = 10 -17 m 2 , ξ = 7×10 -4 and ρ c = 0.39, the mean crack length c can be calculated to be around 0.1 mm, which is in good agreement with direct observations (Figure 1 C and D) and corresponding to the average crystal length. It is now well-accepted that effective pressure has an influence on permeability [e.g. [START_REF] Brace | Permeability of granite under high pressure[END_REF][START_REF] Brace | Permeability from resistivity and pore shape[END_REF][START_REF] Brace | A note on permeability changes in geologic material due to stress[END_REF] and that the effective stress law for the permeability is not straightforward [e.g. [START_REF] Ghabezloo | Effective stress law for the permeability of a limestone[END_REF][START_REF] Braun | Theoretical analysis of pore pressure diffusion in some basic rock mechanics experiments[END_REF]. This is due to the potential closure of pre-existing cracks [e.g. [START_REF] Brace | Permeability of granite under high pressure[END_REF][START_REF] Brace | Permeability from resistivity and pore shape[END_REF][START_REF] Brace | A note on permeability changes in geologic material due to stress[END_REF].

Thus, pressure pulse propagation could lead to variations of permeability due to effective pressure variations.

Combining mass conservation and Darcy's law (equation 3) and neglecting fluid and rock compressibilities, fluid propagation can modelled with a simple 1D diffusion equation [e.g. [START_REF] Rice | Some basic stress diffusion solutions for fluidsaturated elastic porous media with compressible constituents[END_REF]:

∂p(z) ∂z = C(t), ( 10
)
where C is a constant. Neglecting rock compressibility may be a strong assumption, but this can be verified using the fibre optic sensors and comparing their measurements to the pore pressure evolution at specific positions (corresponding to those of the fibre optic sensors) calculated with a constant pore pressure gradient derived from pressure values in the upstream and downstream pumps (Figure 9 Fibre 2 may not be very reliable (section 2.3) and permeabilities calculated with this fibre could not be considered. Permeability calculated between upstream and fibre 1 is lower

D R A F T January 20, 2020, 9:18pm D R A F T X -20 NICOLAS ET AL.: FIBRE OPTIC LOCAL PORE PRESSURE
than permeability calculated between fibre 3 and downstream, implying that permeability could decrease with decreasing effective pressure. Permeability calculated between fibres 1 and 3 (purple points in Figure 9 c) is comprised between the previous permeabilities.

Comparison between these three permeabilities may not be reliable because the fibre optic sensor measurements have a precision that can be estimated to be around 0.5 MPa (see section 2.3). The evolution of permeability as a function of time for one section of the sample is possible to consider. Permeability calculated in the section of the sample nearest to the downstream side increases with pore pressure, due to crack opening with effective pressure decrease due to pore pressure increase (black dots in Figure 9 c). Permeability calculated in the section of the sample nearest to the upstream side remained constant as a function of time (blue dots in Figure 9 c), which can be explained by the fact that pore pressure remained almost constant in this part of the sample. The mean permeability calculated between the upstream and downstream pumps is always between the local permeability measurements. Finally, when pore pressure was almost at equilibrium, local pore pressure measurements were not precise enough to calculate local permeabilities.

However, this pore pressure pulse propagation experiment highlights the interest of the use of fibre optics to measure local pore pressures.

Modelling pore pressure increase in a disconnected dead volume

Following [START_REF] Mayr | Acoustic emission induced by pore-pressure changes in sandstone samples[END_REF], the experiments were modelled considering the diffusion of water from the upstream side, through the sample and into the downstream dead volume

∂p(z, t) ∂t = D ∂ 2 p(z, t) ∂z 2 , ( 11
)
where D is the diffusivity.

Experimental results show that downstream dead volume controls the increase of pore pressure downstream and permeability is relatively high, implying that the gradient is the sample is controlled by the downstream dead volume. Let us examine the upstream and downstream boundary conditions. Initial pore pressure is equal to 0.2 MPa upstream, in the complete sample and in the downstream reservoir:

∀z, p(z = 0, t < t i ) = 0.2MPa, ( 12
)
where t i is the time at which the pore pressure step is applied. Then, even though the upstream pressure application may be slightly delayed, the application of the pore pressure step upstream is simply approximated by a step function, such as:

p(z = 0, t ≥ t i ) = 20MPa. ( 13
)
Let us now examine the role of the downstream dead volume, as done by Pimienta et al.

[2016]. We consider fluid mass continuity downstream, i.e., the change of mass in the sample equals the fluid mass change in the dead volume, with opposite sign. Following [START_REF] Brace | Permeability of granite under high pressure[END_REF] and [START_REF] Lin | Compressible fluid flow through rocks of variable permeability[END_REF], the downstream boundary condition is:

S 1 ∂p d ∂t z=L + kA η ∂p p ∂z z=L = 0, ( 14
)
where L is the sample length, and S 1 is the downstream reservoir storage capacity. The storage capacity is:

S 1 = V 1 K f , (15)
where V 1 is the reservoir volume and K f is the fluid bulk modulus. Note that we neglect the reservoir compressibility. Equations 11 and 14 are resolved with the boundary conditions expressed in equations 12, 13, and 14. Note that this implies that permeability is

D R A F T January 20, 2020, 9:18pm D R A F T X -22 NICOLAS ET AL.: FIBRE OPTIC LOCAL PORE PRESSURE
considered constant during the complete experiment, whatever the effective pressure. We thus neglect the permeability sensitivity to effective pressure.

The values of the parameters used in the model are given in Table 1. Rock porosity, bulk Figure 11 shows the comparison between our experimental results and the model predictions. Predicted downstream pore pressure evolution (black curve) is in very good agreement with the experimental results. However, pore pressure evolution at the three points along the sample corresponding to the location of the fibre optic sensors are in moderate agreement since the present model is not able to predict the sudden increase to the intermediate values followed by a smooth increase to their maximum value. This could be explained by the fact that the applied pore pressure pulse could act as a stress boundary at the bottom of the sample (since it does not diffuse instantaneously). This additional external stress could lead to a compression of the downstream part of the sample.

Due to poro-elastic behaviour, this additional compaction could lead to a sudden pore pressure increase without diffusion, which was not considered by the simulation. Another explanation could be that permeability evolution as a function of effective pressure (and thus pore pressure) has to be taken into account in the model.

Implications and perspectives

The use of fibre optics to measure local pore pressure was shown to be useful to understand the hydro-mechanical behaviour of low-permeability porous materials. In this

D R A F T January 20, 2020, 9:18pm D R A F T NICOLAS ET AL.: FIBRE OPTIC LOCAL PORE PRESSURE X -23
study, fibre optics were used to monitor the migration of a pore pressure pulse but it also opens new horizons for this technique, which are discussed below.

Pore pressure diffusion could also control the transmission of stress and the diffusing pore

Conclusions

It is well known that pore pressure has a major influence on the mechanical behaviour of porous rocks, and thus is of interest for industrial purposes and for our understanding D R A F T January 20, 2020, 9:18pm D R A F T of natural events such as earthquakes. A key parameter for fluid pressure build up is the ability for the material to diffuse fluid, that is to say permeability. However, permeability can be very variable. In the heat-treated andesitic rock used in this study, permeability decreased from 8 × 10 -16 m 2 to 1.5 × 10 -17 m 2 with a confining pressure increase from 2

MPa to 40 MPa. Moreover, a differential stress of 356 MPa decreased permeability from 10 -17 m 2 to 10 -19 m 2 .

In low permeability samples, fluid diffusion can be very slow, implying that it is hard to follow from the sides of the sample, especially when there are large dead volumes. To This calibration data is used to calculate the local pore pressure as a function of the three fibres response during the pore pressure pulse propagation. The inset is a zoom on the three flows when upstream flow is at its minimum value (0.01ml/min).

  meability measurement, following[START_REF] Fortin | Influence of thermal and mechanical cracks on permeability and elastic wave velocities in a basalt from Mt. Etna volcano subjected to elevated pressure[END_REF], we consider that the fluid flow due to the compaction of the sample occurs on both sides of the sample. The rate of fluid volume output due to porosity compaction can simply be calculated as the mean of the fluid flow on both sides of the sample. Upstream imposed flow and downstream flow are significantly higher than compactive volume rate (mean flow) as shown in Figure4, implying that permeability can be calculated with the imposed flux and neglecting the output of fluid due to porosity volume changes.

  Figure 6 (b) shows a linear interpolation of pore pressure along the sample as a function of time. Magnifications of the behaviour at the beginning of the experiment are provided in Figure 6 (c to f)increased upstream, the fibre optic sensors reacted one after the other as a function of their location along the sample (Figure 6). They first showed a sudden increase to intermediate values ranging from 15 MPa (the nearest from upstream) to approximately 5 MPa (nearest to downstream). Downstream pore pressure increase was very smooth (Figure 6 a). Fibre optic pressure measurements remained between the two extremes (downstream and upstream) pore pressures, indicating a pressure gradient along the sample. After 28 minutes, pore pressure was almost constant in the whole sample at values around 19 to 20 MPa.

Figure 7

 7 Figure 7 (a) shows the evolution of differential stress and volumetric strain during tri-

  non-linearity in the stress-strain relation indicates crack closure up to a hydrostatic stress of 28 MPa, leading to a decrease of permeability from 8 × 10 -16 m 2 to 1.5 × 10 -17 m 2 which is similar to previous observations [e.g.[START_REF] Vinciguerra | Relating seismic velocities, thermal cracking and permeability in mt. etna and iceland basalts[END_REF] Guéguen et al., 2011;[START_REF] Nara | Influence of macro-fractures and micro-fractures on permeability and elastic wave velocities in basalt at elevated pressure[END_REF].

5. 2 .

 2 Describing fluid pressure pulse propagation under hydrostatic stress

  a). A zoom into the beginning of the pressure pulse is provided in Figure 9 (b). Calculated pore pressures are close to measured data, showing that the pore pressure evolution follows an evolving constant gradient along the sample during pore pressure increase and decrease, in agreement with equation 10. Assuming a constant flow along the sample (i.e., steady state or neglected fluid compressibility), local pore pressure measurements and Darcy's law were used to calculate the evolution of permeability along the sample as a function of time (Figure 9 c). Local permeability measurements are in good agreement with the mean permeability calculated from upstream and downstream pump pressure measurements (light blue in Figure 9 c).

(

  tubing and pump) as shown in Figure 10. Neglecting fluid compressibility but removing the hypothesis that rock compressibility can be neglected compared to equation 10, fluid propagation can modelled with a new diffusion equation. During pore fluid migration, pore pressure p in the sample depends on the distance from upstream z and time t, which by a 1D version of the diffusion equation [e.g. Rice and Cleary, 1976]:

  modulus and intrinsic permeability are taken equal to those measured on heat-treated samples. Fluid characteristics (bulk modulus, density and dynamic viscosity) are taken equal to those of water commonly obtained at 20 o C. Downstream volume was experimentally measured on the machine. To sum up, the values used for the model are taken equal to those experimentally obtained on rock samples, on the machine and values for fluid (water) are taken from the literature.

  pressure front could play a dual role in the triggering of seismicity by decreasing the coefficient of friction and by decreasing the strength [e.g.[START_REF] Talwani | Pore pressure diffusion and the mechanism of reservoirinduced seismicity[END_REF]. Excess pore pressure can also cause an embrittlement of porous fluid-saturated media, as shown byFarquharson et al. [2016a] in the case of volcanic edifices. Comparing pore pressure migration with fibre optics and acoustic emission migration with an acoustic system could help improving our understanding of triggered seismicity. Rapid fault slip during earthquakes can increase pore pressure because of thermal pressurization, which can decrease shear resistance of the fault material [e.g. Rempel and Rice, 2006; Sulem et al., 2007], or pore pressure decrease because of the opening of new fractures. It would be possible to use fibre optics during failure of samples to have a closer look at pore pressure evolution during shearing inducing dilatancy and failure of low-permeability porous media subjected to triaxial stresses. Fibre optics would be useful to study the evolution of local pore pressure during laboratory stick-slip experiments [e.g. Passelègue et al., 2013] under water saturated conditions, which is similar to what is observed on natural faults [e.g. Passelègue et al., 2014]. Natural faults could control the pathways of fluids and thus fluid diffusion in geothermal fields [e.g. Brehme et al., 2016b]. Fibre optics can be used to monitor the evolution of permeability in low-permeability saw-cut samples to study the evolution of permeability of natural faults as a function of confining pressure and differential stress.

Figure 1 .

 1 Figure 1. Micrographs of intact andesite (A and B) and of 930 o C heat-treated andesite (C

Figure 2 .

 2 Figure 2. Schematic view of the experiment. A pore pressure increase from 0.2 MPa to 20

Figure 3 .

 3 Figure 3.Response of the three fibre optics to an increase of pressure from 2 MPa to 40

Figure 4 .

 4 Figure 4. Imposed upstream flow (blue), downstream flow (red) and mean flow (green) are

Figure 5 .

 5 Figure 5. (a) Evolution of the volumetric strain plotted versus confining pressure during

Figure 6 .

 6 Figure 6. (a) Evolution of controlled upstream (blue) and measured downstream (black) pore

Figure 7 .Figure 8 .

 78 Figure 7. (a) Evolution of differential stress and volumetric strain as a function of time during

Figure 9 .

 9 Figure 9.(a) Evolution of controlled upstream (blue) and measured downstream (black)

Figure 10 .

 10 Figure 10. Schematic view of the model. We consider the diffusion of pore fluid through a

Figure 11 .

 11 Figure 11. Comparison of data and results of the model for pressure measurements upstream,

Table 1 .

 1 Summary of the parameters and their values used in the model.

	Parameter	Symbole	value
	Porosity	φ	1.5%
	Downstream volume	V1	50 ml
	Fluid bulk modulus	K f	2.2 GPa
	Rock bulk modulus	K d	9.5 GPa
	Sample radius	r	2.5 cm
	Sample length	L	10 cm
	Intrisic permeability	k	10 -18 m 2
	Fluid density	ρ	1000 kg/m 3
	Gravitational force	g	9.8 m/s 2
	Fluid dynamic viscosity	η	1.002 × 10 -3 Pa.s
	D R A F T		January 20, 2020, 9:18pm	D R A F T


Maxime Thieffry

Alexandre Kruszewski

Thierry-Marie Guerra

Christian Duriez

Christian Duriez Lpv

LPV Framework for Non-Linear Dynamic Control of Soft Robots using Finite Element Model

Keywords: Soft robotics, Control, Linear parameter-varying systems, Lyapunov methods, Model Reduction

établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

LPV framework for Non-Linear Dynamic Control of Soft Robots using Finite Element Model

Maxime Thieffry * Alexandre Kruszewski * * * Thierry-Marie Guerra * * Christian Duriez * * *

INTRODUCTION

The desire to design robots to interact with humans and/or in confined space in contact with the environment pushed the robotics community to invent new paradigms, among which is soft robotics. Biology provides a major inspiration in the design of soft robots: the elephant trunk, octopus arm or snake body lead to the design of the first soft robots [START_REF] Majidi | Soft robotics: a perspective -current trends and prospects for the future[END_REF][START_REF] Kim | Soft robotics: a bioinspired evolution in robotics[END_REF]. By soft robots, we mean robots made of deformable materials whose motion are obtained by deformation and where no joints are presents in the structure. Soft robots have an theoretical infinite number of degrees of freedom. Developing a mathematical model suitable to describe the dynamics of such systems is a challenging task as it must be at the same time computationally affordable and sufficiently accurate.

Different methods have been proposed to solve this challenge, like constant curvature [START_REF] Webster | Design and kinematic modeling of constant curvature continuum robots: A review[END_REF]. This method eases the modeling of soft structure due to geometrical assumptions that restrict the design of the robots to beam-like structures. It has been successfully applied to many continuum robots with different actuation systems such as cables or pneumatic chambers [START_REF] Marchese | Design and control of a soft and continuously deformable 2d robotic manipulation system[END_REF]. Initially developed for kinematics studies, this method has been extended to dynamics modeling and used in closed-loop control experiments in [START_REF] Katzschmann | Dynamic motion control of multi-segment soft robots using piecewise constant curvature matched with an augmented rigid body model[END_REF].

In addition, Cosserat theory provides a geometrically exact method to model soft robots, taking into account large deformations and displacements and handling material nonlinearities [START_REF] Trivedi | Geometrically exact models for soft robotic manipulators[END_REF]. As for the constant curvature method, it is intended for beam-like robots.

Numerical methods have been proposed to model soft robots without limitations on the geometry. The Finite Element Method (FEM) consists in discretizing the structure into a finite number of small elements, the underlying equations coming from continuum mechanics are then solved for each of these elements. A modeling and simulation software dedicated to soft robots is implemented upon the open-source framework SOFA along with a plugin for the specific needs of soft robots [START_REF] Coevoet | Software toolkit for modeling, simulation, and control of soft robots[END_REF]. This modeling method can handle most of the geometries, provided that a CAD file for the structure exists.

To design feedback controllers for soft robot, depending if one is interested in controlling the kinematics or dynamics of the structure, it can lead to different controller designs. Many different methods exist to tune a closed-loop controller, depending on the modeling strategy, the targeted application or an optimization criterion. If some results have shown significant results for kinematics control, few methods exist for dynamic control (see [START_REF] Thuruthel | Control strategies for soft robotic manipulators: A survey[END_REF]).

Kinematics controllers may not be sufficient to perform high speed tasks (such as jumping) or high speed obstacle avoidance. Recently, different authors have proposed new methods to control the dynamic behavior of soft robots. A dynamic controller based on constant curvature model is presented in [START_REF] Falkenhahn | Model-based feedforward position control of constant curvature continuum robots using feedback linearization[END_REF] where the model is used to generate a feed-forward action coupled with a PID controller. Also based on piece-wise constant curvature model, authors of [START_REF] Della Santina | Dynamic control of soft robots interacting with the environment[END_REF] present a dynamic controller that enables dynamic trajec-tory tracking for a continuous soft robot while handling interactions with the environment. Combining Cosserat and Lagrange dynamic models with Ritz-Galerkin methods, [START_REF] Sadati | Control Space Reduction and Real-Time Accurate Modeling of Continuum Manipulators Using Ritz and Ritz-Galerkin Methods[END_REF] presents a real-time model of continuum manipulators that enables nonlinear impedance and configuration control. In addition, model predictive control has been used recently to perform trajectory following tasks [START_REF] Bruder | Modeling and control of soft robots using the koopman operator and model predictive control[END_REF]. Authors use Koopman operator and system identification methods to construct a discrete-time linear model based on which a linear model predictive controller is designed.

CONTRIBUTIONS

This paper aims at extending our previous work by considering the non-linear model of the robot for the control design. In [START_REF] Thieffry | Control design for soft robots based on reducedorder model[END_REF], we have presented a control strategy based on a finite element model of soft robot dynamics. and it has been validated on experimental platform in (Katzschmann et al., 2019a). However, the control law was designed using a linear time invariant (LTI) model, this is a known limitation of this work as it constrains the guaranteed domain of the controller. This papers presents a methodology to build reduced order linear parameter varying (LPV) models of soft robots based on a finite element model. A linear controller is then designed for this LPV model.

FINITE ELEMENT MODEL

Linear large-scale model

To overcome the difficulty of providing realistic analytical models of soft robots, we rely on numerical methods to solve the equations coming from continuum mechanics. The finite element method is commonly used in computational mechanics and is implemented within the SOFA framework [START_REF] Coevoet | Software toolkit for modeling, simulation, and control of soft robots[END_REF] to obtain real-time performances. This modeling method consists of discretizing the deformable solid into small elements. The equations are then solved for each of these elements and assembled to obtain the behavior of the hole structure of the solid.

Let us start with the formulation given by the second law of Newton that models the dynamic behavior of a body as: M(q)q = P(q) -F(q, q) + H T (q)λ (1) where q ∈ R n is the vector of generalized degrees of freedom, M(q) : R n → R n×n is the inertia matrix, F(q, q) : R n × R n → R n represents the internal forces applied to the structure and P(q) : R n → R n gathers known external forces. Finally, matrix H(q) : R m → R n×m is the matrix containing the actuation directions while λ ∈ R m is the vector of actuators forces.

Let q 0 ∈ R n be a stable equilibrium point. It is induced by the gravity field P(q 0 ) and the actuation input λ 0 . As q 0 is an equilibrium point, it holds: 0 = P -F(q 0 , 0) + H T (q o )λ 0

(2)

The internal force vector F is a nonlinear function of the positions and the velocities. We apply a Taylor series expansion to F and make the following first order approximation around this configuration q 0 and it holds:

F(q, v) = F(q 0 , 0) - ∂F(q, v) ∂q q=q0 v=0 ∂q - ∂F(q, v) ∂v q=q0 v=0 ∂v (3)
We define the stiffness matrix K and the damping matrix D as:

K : R n × R n → R n×n K(q, v) = ∂F(q, v) ∂q D : R n × R n → R n×n D(q, v) = ∂F(q, v) ∂v (4)
Let us define the vector of velocity v ∈ R n , v = q, the state x ∈ R 2n of the robot is then defined as:

x = v q (5)
In the following, we assume that the external forces P(q), the mass matrix M(q), the stiffness K(q, v), the damping D(q, v) and the actuation direction H(q) are constant around q 0 :

P = P(q) = P(q 0 ) H = H(q) = H(q 0 ) M = M(q) = M(q 0 ) K = K(q, v) = K(q 0 , 0) D = D(q, v) = D(q 0 , 0) (6)
Finally, we define the displacement vector d as d = q -q 0 (7)

From ( 1) and ( 3), the equation of motion around this equilibrium point q 0 writes:

M v = P -F(q 0 , 0) -Kd -Dv + H T λ (8)
Computing ( 8)-(2) yields to the following equation, modeling the motion around an equilibrium point q 0 :

M v = -Kd -Dv + H T (λ -λ 0 ) (9)
Under this assumption and defining a new input u = λλ 0 , we obtain a linear model of the behaviour of the robot described by the following equation:

M v = -Kd -Dv + H T u (10)
Considering the state vector x = (v d) T and provided that the matrix M is regular, equation ( 10) can be written as a Linear Time Invariant (LTI) model:

   ẋ = -M -1 D -M -1 K I 0 x + -M -1 H T 0 u y = Cx ( 11
)
where I is the identity matrix of dimension n and the matrix C is a sparse matrix defining the end-effector coordinates.

Reduced Order Model

In equation (1), the vector q is the displacement of each nodes of the mesh in the three dimensions of space x, y and z. The more precise the model is, the more variables there are in the model. For soft robotics applications, the mesh of the robot is often made of thousands of variables, that makes system (11) a model of large dimensions. To overcome the dimensionality issue arising when studying these systems, model order reduction is used in this work. Consider a nonlinear model: 12) projection-based model order reduction consists of decomposing the full-order state x into two parts, a low-order state x r ∈ R r and a neglected state x r ∈ R n-r such that:

ẋ(t) = f (x(t), u(t)) , x ∈ R n , n 1 (
x = V r x r + V r x r with x r = W T r x x r = W T r x ( 13
)
The problem is thus to find two projectors V r ∈ R n×r and W r ∈ R n×r to compute a x r with r n such that:

x r = W T r x ; x ≈ V r x r (14)
In other words, the approximation method consists of finding the matrix V r and W r such that:

ẋr (t) -W T r f (V r x r (t), u(t)) ≈ 0 (15)
This provides an approximation of the large-scale system:

ẋr (t) ≈ W T r f (V r x r (t), u(t)) , x r ∈ R r , r n (16)
Singular Value Decomposition (SVD) based methods and Krylov (moment-matching) based methods are commonly used in the literature [START_REF] Benner | Model Reduction and Approximation: Theory and Algorithms[END_REF]. Proper Orthogonal Decomposition (POD) is a SVD-based method that is directly usable for non-linear large-scale systems, that makes is suitable for our application.

The underlying idea behind POD algorithm is to collect samples of the state of the studied system and find a reduced basis that approximates these samples. This method is well suited for applications where simulation is available, as it is easy to get a collection of samples of states. The SOFA framework used to model soft robots includes algorithms to perform model reduction using POD [START_REF] Goury | Fast, generic, and reliable control and simulation of soft robots using model order reduction[END_REF].

Let Σ : ẋ = f (x(t), u(t)) be a nonlinear system and let S be a collection of s snapshots.

S = (x(t 1 ) x(t 2 ) . . . x(t N )) ∈ R n×s
Computing a SVD of this matrix, we get V , the left singular matrix of S :

V ΣΩ T = S (17)
Let us define V r as the r th first columns of V , we have the following approximation: Then, the projectors V r and W r are simply obtained as:

S = V ΣΩ T = V r Σ r Ω T r + ∆ (
V r = W r = V r (19) so that it holds:

x r = W T r x ; x ≈ V r x r (20) In addition, projectors are orthogonal and it holds:

W T r V r = I r and W T r V r = 0 (21)
Applying POD to the linear model ( 11), the reduced order model writes:

ẋr = W T r AV r x r + W T r Bu = A r x r + B r u y r = W T r Cx r = C r x r (22)
4. ILLUSTRATIONS

Trunk-like robot

This experimental platform is a trunk-like robot presented in figure 1 with a schematic view in figure 2. It is made entirely of silicone, it is 18 centimeters long and the thickness at its base and its tip are respectively 2.5 and 1 centimeters. The structure is driven by 4 cables -actuated by 4 servomotors whose entry are the cable length -that permit to work in the 3 dimensions of space. The output of the system is the position of the tip (red point in figure 1) in the three dimensions of space x, y and z. A comparison of different FEM mesh with different accuracy is given in figure 3. A mesh with 210 nodes is not accurate enough to represent the geometry faithfully.

Conversely, a mesh with 6012 nodes does not give more accurate results compared to a medium size mesh of 1557 elements. This medium size mesh is a good compromise and will be used in this work. With 1557 nodes in the mesh, the dimensions of the state vector in system ( 11) is also 1557 × 3 × 2 = 9342 state variables (3 directions of space for displacement and velocity). The output is the position of the end-effector (red-point of figure 1) in the three directions of space.

A linear controller has been designed in [START_REF] Thieffry | Reduced Order Control of Soft Robots with Guaranteed Stability[END_REF] where results show how to control the dynamics of the robot in simulation experiments. Simulation and realworld experiment shows the effectiveness of the approach. The algorithm proposed in previous work is able to control the dynamic behavior of the soft robots. However, as the model was linear, the workspace where the algorithm has been validated was limited to a region around the equilibrium point. In addition, when the target position is outside of the region of validity of the linearization assumption, the results are unpredictable and can lead to undesired behavior, as shown in figure 4. This limitation may have several explanations but it is probably due to the fact that the target position is outside the validity range of the linearization. Next section presents tracks to model this kind of soft robots using LPV models.

LPV MODEL OF SOFT ROBOTS

Basic Concepts of LPV models

Linear Parameter varying (LPV) systems have gained popularity during the 1990s, benefiting from the extension of H ∞ optimal control. Many nonlinear systems can be written as quasi-LPV systems, which has two advantages: the first one is to avoid writing a nonlinear model that requires precise knowledge of the process studied. The second one is to take advantage of all the techniques developed for LTI systems. Polytopic systems are a common way of modeling LPV systems [START_REF] Apkarian | Parameterized lmis in control theory[END_REF] and are of the form:

A(ρ(t)) B(ρ(t)) C(ρ(t)) D(ρ(t)) = N i=1 ρ i (t) A i B i C i D i (23)
The birth of LPV systems comes from gain scheduling techniques, where the idea is to linearize nonlinear systems around different operating points yielding to a collection of local LTI models. Then, interpolation functions link each

Multiple linearizations along a trajectory

To extend the guaranteed domain of the control algorithms, a solution is to design a controller valid for the nonlinear system (1), i.e. without assumption about area of validity. With x = (v, q) T , this model writes:

ẋ = A(x)x + B(x)u (24)
This defines a non-linear large-scale state-space system. While it is already challenging to design a controller for nonlinear systems, the complexity of the problem increases with the dimensions of the system.

Let us consider a fixed number K of equilibrium point x e k , k ∈ {1, . . . , K} along the trajectory defined in figure 5. These equilibrium points are induced by the gravity field and a collection of inputs u e k , such that:

0 = A(x e k )x e k + B(x e k )u e k ( 25
)
Around each of the equilibrium point x e k , let us define z k = x -x e k and u k = u -u e k to write the collection of linear systems:

żk = A k z k + B k u k , k ∈ {1, . . . , K} (26)
The POD reduction algorithm is well suited for nonlinear systems. The snapshots are captured so that the entire workspace of the robot is modeled, in other words all the subsystems (A k , B k ) are included in the snapshot space.

The reduction provides the projection matrices V r and W r that are valid for all of the subsystems. For the entire workspace of the robot, the projection writes:

x r = W r x ; x ≈ V r x r (27)
The key-point of using the POD algorithm, is that all subsystems (A k , B k ) share the same projectors, yielding to a collection of low order linear systems:

żr k = A r k z r k + B r k u i ( 28
) where z r k = x r -x re k and with k ∈ {1, . . . , N }.

Approximation of the collection of linear systems

Let us study the evolution of the coefficient of the reduced order system matrix A r of the cable-driven soft robot (see figure 1). The reduction is done using POD algorithm and the reduced order system is of dimension 6. The matrix A r also writes:

A r k =    a 11 k a 12 k . . . a 16 k a 21 k . . . . . . a 26 k a 61 k . . . a 65 k a 66 k    (29)
The coefficient of the matrix A r are saved at 23 different locations (i.e. K = 23) along the trajectory depicted in figure 5. The evolution of some of the matrix coefficients is shown in figure 6. One interesting track is to use a sparse approximation technique to reconstruct all of the coefficients from matrices A r and B r with a minimum of variables. According to the nice properties of radial basis functions, such as sparse universal approximation, these first trials are focused on these functions.

A radial basis function is a real-valued function f whose value depends only on the distance from the origin, so that f (x) = f ( x ); or alternatively on the distance from a given point c, so that f (x, c) = f ( x -c ). This point c is then called a center. Gaussian function of the form f (x) = e -x 2 are commonly used.

A radial basis network is created to approximate a function defined by a set of data points, in our case the coefficients (entries) of each matrices coming from the linearized systems along the trajectory. A radial basis network is a network with two layers, a hidden layer of radial basis neurons and an output layer of linear neurons. The MATLAB functin newrb is used to create the radial basis functions network.

The inputs of the network are the coefficient of the systems matrices (A r k , B r k ) for each subsystems, the networks output is the function that approximates those coefficient that writes:

λ i (φ) = e -(bi φ-ci ) 2 , i ∈ {1, . . . , N } (30
) where N is the number of functions λ used to approximate the matrix coefficients. In addition, b i are the bias values, φ are the input weight values and c i are the centers values.

Construction of LPV model

From this approximated function, we get the following system:

ẋr = N i=1 λ i (A ri x r + B ri )u (31)
However, due to the construction of the radial basis function network, only λ > 0 holds and there is no reason for the functions to share the convex sum property, i.e.

N i=1 λ i = 1. Nevertheless, the sector nonlinearity approach is a systematic way to derive a polytopic model [START_REF] Tanaka | Fuzzy control systems design and analysis: a linear matrix inequality approach[END_REF]. Values of λ i are bounded and it holds λ i ∈ [m, m] ∈ [0, 1], i ∈ {1, . . . , N } , one can also write λ i = ω i 0 m + ω i 1 m, where:

ω i 0 = m -λ i m -m ; ω i 1 = λ i -m m -m (32)
and the function ω i i satisfy the convex sum property: ω i 0 + ω i 1 = 1 and ω i i ≥ 0. Let us define functions h i (Bernal and Guerra, 2010):

h 1+i1+i2×2+...+i N ×2 N -1 = N j=1 ω j ij (33)
We finally end with a polytopic model ( 34), based on which the design of a controller is detailed in next section. ẋr =

2 N i=1 h i ( Ãri x r + Bri u) (34) where Ã1+i1+i2×2+...+i N ×2 N -1 = λ 1 (φ i1 )A r1 + • • • + λ N (φ i N )A r N = j=1...N (λ j (φ ij )A rj ) (35)
with φ 0 = min(φ) and φ 1 = max(φ).

CONTROL DESIGN

Theory

The simplest control design is a linear state feedback controller u = -Lx r (36) corresponding to the closed-loop system:

ẋr = 2 N i=1 h i ( Ãri -Bri L)x r (37)
The closed-loop LPV polytopic system (37) is quadratically stable if and only if there exists a matrix P = P T > 0 such that: ( Ãri -Bri L) T P + P ( Ãri -Bri ) < 0 (38) hold for all i = 1, . . . , 2 N , see [START_REF] Boyd | Linear matrix inequalities in system and control theory[END_REF].

Condition (38) corresponds to a set of N LMI constraints that, depending on the number of vertices of the polytope, can sometimes lead to a huge number of constraints to solve. As the number of parameters N growths, so does the complexity of the condition to satisfy. It can also be time and memory consuming as the number of LMIs to solve is 2 N .

APPLICATION TO STUDIED SOFT ROBOT

Approximation using RBF

The evolution of the coefficients along the trajectory studied in this work is done using the RBF network that Then, a second experiment aims at driving the robot from its rest shape to a highly deformed position, where a linear model is not able to capture the dynamic behavior. The robots starts at its rest shape (figure 5 a) and converges to a highly deformed position (figure 5 c). Results are shown on figure 9. Compared to results of previous works shown in figure 4, figure 9 shows that designing a controller based on LPV models makes it possible to control the robot in a wider workspace.

DISCUSSION AND FUTURE WORK

Simulation experiments shows the validity of the approach using a linear controller for the LPV system. Better performances could be reached with more elaborate control laws, such as the parallel distributed compensation (PDC) [START_REF] Tanaka | Fuzzy control systems design and analysis: a linear matrix inequality approach[END_REF]. This is the topic of ongoing research. It is composed of linear feedbacks assembled together with the same nonlinear function h i as in the model ( 34):

u = - 2 N i=1 h i L i x (39)
yielding to corresponding closed-loop model: In addition, to study the robustness of the control law, one should study the modeling errors coming from the reduction algorithm. Indeed, the POD algorithm decompose the full order state x into two parts, a reduced part x r that represents faithfully x and a neglected part x r that is neglected such that: x = V r x r + V r x r ≈ V r x r (41) Results of [START_REF] Thieffry | Control design for soft robots based on reducedorder model[END_REF] shows how to design a controller that takes into account this reduction error x r based on linear models. Its extension to LPV models has to consider a global error coming from the convex aggregation of the polytopic vertices. This extension is left for future research.

CONCLUSION

This work presents new results about dynamic control of soft robots using a nonlinear reduced-order numerical model. The finite element method provides us with a largescale nonlinear model that is linearized along a trajectory to build a linear parameter varying (LPV) model. This work is an important step towards nonlinear control of soft robots modeled using finite element method. Further

  18)where Σ r contains the r-first singular values of S and ∆ represents the model order reduction errors.

Fig. 1 .

 1 Fig. 1. Soft robot used for experimental validation, fully made of silicone. Top: front view of the robot. Bottom: side view. Red = end-effector, location of the sensor.

Fig. 2 .

 2 Fig. 2. Design of the robot: slice view (left) and side view (right).The robot is actuated with 4 cables in red.

Fig. 3 .

 3 Fig. 3. Comparison of different FEM mesh of the Trunklike robot with different accuracy. Top: FEM mesh with 210 nodes, Middle: FEM mesh with 1557 nodes, Bottom: mesh with 6012 nodes.

Fig. 4 .

 4 Fig. 4. Unsuccessful attempt to reach a more target distant from the linearization point. From top to bottom: displacement of end-effector along x, y and z axis in red with the desired value in black.

Fig. 5 .

 5 Fig. 5. Trajectory along which the model is linearized. local subsystems. These interpolation functions are called scheduling function and describe the change of operating point. The scheduling signal is referred to as ρ. Therefore, the resulting controller are dependent on the varying signal ρ, thus resulting in parameter varying systems.

Fig. 6 .

 6 Fig. 6. From left to right, top to bottom : evolution of coefficient a 11 , a 14 , a 21 and a 25 of matrix A r along the trajectory, from equation (29).

Fig. 7 .

 7 Fig. 7. Comparison of matrix coefficients a 11 , a 14 , a 21 and a 25 along the trajectory (red dots) and the outputs of the radial basis network for these coefficients (blues lines).provides 4 functions λ, i.e. N = 4. The mean square error (MSE) between the network input (matrix coefficients along the trajectory) and the network output for the studied soft robot is M SE = 0.39 given that the maximum values for a coefficient along the trajectory is 7.85. This represents a maximum error of 4.97%. A comparison of the matrix coefficients a 11 , a 14 , a 21 and a 25 and the output of the radial basis network for these coefficients is shown in figure7.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Displacement of end-effector along x (top) y (center) and z axis (bottom) in open-loop (dashed line) and closed-loop (plain line). The robot starts from its rest shape (Fig. 5 a) and converges to a slightly deformed position (Fig. 5 b).

This research was conducted with support of ANR (Project ANR-17-ERC2-0029), the European Union through the European Regional Development Fund (ERDF), the French Ministry of Higher Education and Research, the National Center for Scientific Research (CNRS), and the Hauts-de-France Region.