© CFMS-CFGI-CFMR-CFG, Published by EDP Sciences 2016

1 Introduction

La recherche sur les circulations de fluides dans les milieux poreux et fracturés a connu ces dernières décennies un essor croissant motivé par des enjeux de sécurité et d’efficacité liés à la réalisation d’ouvrages en souterrain (mines et carrières, stockages) et à la surface du sol (fronts rocheux, barrages). Ces ouvrages sont notamment liés à la problématique de l’enfouissement des déchets (nucléaires, industriels et ménagers), de la production ou du stockage d’énergie, ou encore de la lutte contre les effets de serre (stockage du CO₂).

La stabilité mécanique et le régime d’écoulement des eaux souterraines sont affectés par le comportement hydraulique des terrains. Or, ce dernier est régi par l’état de fracturation du massif rocheux, mais aussi par la porosité et la perméabilité de la matrice rocheuse.

On conçoit l’intérêt de prendre en compte ces deux aspects, aussi bien lors d’expérimentations en laboratoire ou in situ que dans une modélisation. Des tentatives existent notamment dans le domaine du génie pétrolier, mais le concept de massif rocheux à double porosité a rarement été employé en géotechnique appliquée aux ouvrages souterrains.

Dans le cadre de cet article, on considérera le milieu rocheux comme la superposition d’un réseau de fractures dite « principales » (ou « clés ») et d’un « milieu poreux continu équivalent ». Ce milieu « équivalent » doit rendre compte indirectement du comportement hydromécanique des fractures secondaires et de la matrice. Les caractéristiques hydromécaniques du milieu équivalent doivent être évaluées par des techniques dites d’homogénéisation (ou d’« upscaling »). Cette approche est motivée notamment par la difficulté de prendre en compte explicitement dans un modèle l’ensemble des hétérogénéités existantes dans le milieu naturel (Blum et al., 2009 ; Pouya et Fouché, 2009).

Ces caractéristiques équivalentes sont déterminées selon une démarche d’homogénéisation numérique mise au point antérieurement (Ababou, 1991 ; Min et al., 2004 ; Oda, 1986 ; Thoraval et Renaud, 2004). Il s’agit tout d’abord d’identifier la géométrie de chacun des « blocs » à homogénéiser (matrice découpée par des fractures secondaires). On calcule ensuite, par simulations numériques, la réponse de ces « blocs » à différents chargements hydromécaniques, de sorte à déterminer les tenseurs de perméabilité et de rigidité d’un milieu poreux continu « équivalent » au milieu fracturé considéré.

La mise en œuvre de ce type de modèle nécessite de répondre préalablement aux questions suivantes : quelles sont les fractures « clés » à considérer pour un site donné (du point de vue de la prédiction des écoulements ou des mouvements de terrain) ? Quel écart observe-t-on dans les prédictions suivant que les fractures supplémentaires sont introduites explicitement, ou à travers une représentation de type « milieu poreux continu équivalent » ? La prise en compte du comportement poromécanique de la matrice rocheuse affecte-t-il significativement la réponse du modèle, ou est-il possible de le négliger ?

Dans la présente étude, nous avons examiné l’impact de diverses simplifications sur la prédiction du comportement hydromécanique d’un massif rocheux fracturé et poreux soumis à des sollicitations hydrauliques. Une approche numérique tridimensionnelle a été développée prenant en compte à la fois le comportement hydromécanique de la matrice rocheuse poreuse et des fractures. Le résultat des simulations a été comparé avec les mesures expérimentales réalisées lors de tests d’injection et de pompage de longue durée effectués dans le site de Coaraze (Cappa, 2006 ; Cappa et al., 2005 ; Guglielmi, 1999). L’analyse des différentes options de modélisation réalisées nous conduit à préciser certaines recommandations méthodologiques.

2 Choix d’un jeu de données pertinent pour l’étude méthodologique

2.1 Présentation du site expérimental

L’étude que nous présentons ici est basée sur des expérimentations réalisées sur le site localisé à proximité du village de Coaraze à 30 km au nord de Nice. Il est situé en rive gauche d’un petit vallon orienté N40. Le vallon marque la limite entre les calcaires du Jurassique en rive droite et les faciès du Crétacé, en rive gauche (Fig. 1). Trois faciès litho-stratigraphiques peuvent être distingués dans le Crétacé (Guglielmi, 1999) :

• –

  un faciès calcaire basal Bérriasien-Valanginien, séparé des calcaires jurassiques par un niveau argilo-sableux glauconieux de 1,5 m d’épaisseur ;

• –

  un faciès marno-calcaire Hauterivien-Barrémien à nouveau séparé du précédent par un niveau argilo-sableux glauconieux de 2 m d’épaisseur ;

• –

  un faciès marneux avec de fines intercalations calcaires d’âge cénomanien.

La partie nord-est du site correspond à une structure monoclinale, réglée par la direction N40-45° et le pendage 45–50°SE des couches. Cette zone est recoupée par une faille inverse N50° à N80°, 70°N dont le rejet correspond à un abaissement relatif du compartiment sud par rapport au compartiment nord. Dans la partie sud-est du site, les couches sont intensément découpées par des failles de direction moyenne N140° et de pendage subvertical à 80°NE. Ces failles délimitent des lentilles tectoniques de 10–15 m de long qui correspondent à des compartiments de plus en plus abaissés vers l’aval du vallon. À l’extrême sud-ouest de la zone, le vallon recoupe un contact majeur du Jurassique sur le Cénomanien.

Ce site a été choisi pour l’étude in situ des processus hydromécaniques dans des roches carbonatées fracturées (Guglielmi, 1999). Il correspond à un petit réservoir fracturé (30 m × 30 m) dans un affleurement de surface qui est constituée d’une succession de bancs calcaire à faible porosité sur 15 m d’épaisseur. Les limites hydrauliques du réservoir sont bien définies, avec des couches imperméables à sa base et son sommet (couche d’argile sableuse – voir Fig. 1), une faille imperméable située en aval vers le sud-ouest (de direction N140°). En amont, vers le Nord, une faille perméable permet d’alimenter en continu le réservoir qui est, par ailleurs, drainé par une source située à sa base. Cette source a été équipée d’une vanne qui a permet le contrôle du niveau piézométrique dans le réservoir (Fig. 1).

Lors de la mise en charge, la quasi-absence de réaction d’un piézomètre situé dans le banc glauconieux basal prouve l’imperméabilité de ce niveau. Quand le niveau piézométrique monte, une source temporaire apparaît au toit des calcaires sous le contact avec le deuxième niveau glauconieux sommital. Ce deuxième niveau joue donc aussi un rôle d’imperméable. L’aquifère est donc bien délimité dans les calcaires perméables bérriasiens-valangiens situés entre les deux niveaux imperméables glauconieux. Pour provoquer une mise en charge suffisante dans ces calcaires, toutes les petites fissures du versant situées en rive gauche du vallon ont été colmatées par du ciment. Ce versant joue donc le rôle d’une limite imperméable sur une dizaine de mètres de hauteur au-dessus de la source.

Fig. 1 Schématisation du site expérimental de Coaraze (Guglielmi, 1999). Description of Coaraze experimental site (Guglielmi, 1999).

2.2 Description de la fracturation du site

Les études préliminaires menées sur le site par Guglielmi (1999) ont montré que le réseau de fractures du site de Coaraze se compose de trois familles principales de fractures : trois failles, espacées de 2 à 3 m, ayant un azimut N50 à N70 et un pendage de 70 à 90° vers le NW ; onze failles ayant un azimut N120 à N140, un pendage de 75 à 90° vers le NE et environ 2 m d’espacement ; douze plans de stratification ayant un azimut N40, un pendage de 45° vers le SE et un espacement de 0,5 à 1,0 m.

Des relevés complémentaires ont été réalisés à cette fin, à l’aide d’un tachéomètre haute précision. Le principe de la mesure consiste à positionner le tachéomètre face au versant rocheux étudié et à mesurer les coordonnées X, Y et Z de points positionnés sur les traces visibles des fractures sur le versant. Près de 500 mesures (visées) ont été réalisées, sachant que certaines fractures ont été caractérisées par plus de 20 points. Chaque campagne de mesure à permis l’amélioration progressive du modèle géométrique, les nouvelles mesures étant souvent conditionnées par la confrontation périodique du modèle avec les observations de terrain. Le résultat final est une cartographie tridimensionnelle exhaustive des principales fractures découpant le réservoir.

2.3 Description des essais réalisés et des mesures hydromécaniques retenues pour l’étude

Le premier type d’expérimentation a consisté à imposer une sollicitation globale au massif par fermeture de la vanne installée en pied de versant (Fig. 1). Quand la vanne est fermée, l’eau monte assez rapidement dans les fractures du massif sur une hauteur d’environ 10 m (du fait de l’imperméabilisation de la partie basse du versant). L’ouverture et la fermeture de la vanne permettent de contrôler le niveau de pression de l’eau dans les fractures. Il apparaît notamment que la pression d’eau se stabilise en moins de 30 minutes dans les fractures les plus conductrices (cette stabilisation est plus lente dans les blocs, moins perméables, découpés par ces fractures). Les résultats de cette expérimentation et leur interprétation sont décrits dans (Cappa et al., 2005). Des compléments d’interprétation sont apportés par Guglielmi et al. (2008) où il apparaît que la réponse du modèle global est très sensible à la prise en compte explicite des discontinuités subverticales présentes sur le site. L’homogénéisation des autres fractures (y compris des plans de stratification) par un milieu continu équivalent est largement acceptable, les fractures verticales expliquant 90 % des mesures réalisées. Toutefois, cette représentation peut conduire à des erreurs importantes localement.

Le deuxième type consiste à imposer une sollicitation plus locale en réalisant un essai d’injection dans une faille subverticale F₂ (Fig. 2). On impose une pression ou un débit dans une fracture et on mesure les variations des pressions hydrauliques et des déplacements normaux au plan de fracture, induites par la sollicitation. Cela a nécessité la réalisation préalable de deux forages horizontaux, espacés d’environ 1 m, recoupant la fracture choisie. On distingue des injections de courte durée (« pulse » quelques secondes) et des essais de pompage et d’injection de plus longue durée (succession de paliers de 30 minutes à une heure). Les injections et les pompages d’eau sont réalisés avec une pompe à membrane connectée par un train de tiges à la chambre. Cette pompe peut imposer des variations de pression allant de 0 à 0,7 MPa. Toutes les expérimentations sont réalisées en condition hydraulique saturée (correspondant à l’état d’équilibre du site à vanne fermée). En ce qui concerne les expérimentations

Dans chaque forage, la faille a été isolée avec un packer gonflable (Fig. 3), afin de créer une section isolée de 0,4 m de long. Le dispositif de mesure est composé des capteurs à fibre optique mesurant simultanément la pression de l’eau et le déplacement normal de la faille. Ce dispositif a été adapté à partir de l’appareil BOF-EX développé par RocTest-Telemacs (Cappa et al., 2006). Le dispositif de mesure est capable de réaliser des mesures à haute fréquence [120 Hz] et avec une grande précision (±1 kPa pour la pression et ± 10⁻⁷ m pour le déplacement).

Nous avons rassemblé (Tab. 1) les mesures obtenues dans la fracture en HM₁ et HM₂ au cours des essais de pompage et d’injection par paliers de longue durée qui seront utilisée dans cet article. Les valeurs imposées sont les valeurs en gras, soit la première colonne (débit) dans le cas des essais de pompage et la deuxième colonne (pression) dans le cas des essais d’injection. Les autres valeurs, mesurées in situ, sont induites par la sollicitation.

La pression hydraulique indiquée est une variation de pression par rapport à un état initial de 43 kPa en HM₁ et de 33 kPa en HM₂ :

• –

  il faut imposer un débit de pompage de 0,77 l/s pour que la pression s’annule en HM₁ (elle vaut alors 22 kPa en HM₂). Le pompage induit des déplacements normaux (fermeture) maximaux de 0,9 mm en HM₁ et de 0,2 mm en HM₂ ;

• –

  les paliers de pressions imposées au cours de l’essai d’injection induisent des débits injectés respectifs d’environ 0,7 l/s et 1,2 l/s en HM₁ et des déplacements normaux (ouverture) respectifs de 1 μm et 2,35 μm. Seule la variation de pression induite a pu être mesurée en HM₂, le capteur de déplacement étant temporairement hors d’usage.

Fig. 2 Description des expérimentations d’injection et de pompage au sein de la fracture verticale F2. Description of the experimental device used for local injections or pumping within fracture F2.

Fig. 3 Sonde de mesure mise en place dans les forages pour mesurer la pression hydraulique et les déplacements normaux dans la faille F2. Experimental device inserted into the boreholes to measure hydraulic pressure and normal displacements within the fault F2.

3 Construction de modèles numériques pour représenter les essais d’injection/pompage par paliers de longue durée

La démarche à consisté à construire, dans un premier temps, un modèle géométrique de petite taille (cube de 8 m de coté dont la projection est représenté sur la Fig. 2) autour du point d’injection (HM₁) à partir des mesures de fracturation faites sur le terrain. Les fractures introduites dans ce modèle géométrique sont :

• –

  la fracture instrumentée F₂ qui est sollicitée (elle intersecte les deux sondages HM₁ et HM₂ respectivement à 4,4 m et 4,6 m du versant) ;

• –

  les joints de stratifications S₁₀, S₉, S₈ et S₇ ;

• –

  les diaclases R₂₀ (localisée en dessous de S₉ et limitée par celle-ci), R_9b (localisée au-dessus de S₉ et limitée par celle-ci) et R_9a (non limitée), R_10b (localisée au-dessus de S₉ et limitée par celle-ci) ; J₁ (comprise entre S₆ et S₇), J₂ (comprise entre S₆ et S₇) et D₆ (non limitée).

Un modèle hydromécanique a ensuite été conçu à partir de ce modèle géométrique. Nous avons tout d’abord privilégié une approche discrète en faisant l’hypothèse que les écoulements se produisaient uniquement à travers les fractures prises en compte dans le modèle. Nous avons ensuite développé une approche mixte, où seules les fractures principales sont prises en compte explicitement, les fractures secondaires et la matrice poreuse étant représentées par un milieu poreux équivalent. Cela a nécessité l’utilisation d’une méthode d’homogénéisation pour calculer les tenseurs de perméabilité et de rigidité de ce milieu poreux équivalent. Le nouveau modèle a permis d’évaluer l’impact de la porosité matricielle et des fractures secondaires sur la simulation des expérimentations locales.

3.1 Modèle hydromécanique 3DEC ne prenant en compte que la porosité des fractures

Peut-on négliger d’introduire dans le modèle les fractures autour de la fracture instrumentée ? Quelles erreurs de prédiction commet-on si l’on n’introduit pas les joints de stratification (moins perméables que les autres fractures), si l’on n’introduit aucune autre fracture que la fracture instrumentée ? Pour répondre à ces questions, trois variantes ont été considérées (Fig. 4). Ces variantes se distinguent par le nombre de fractures introduites dans le modèle : (géométrie A) uniquement la fracture instrumentée (F₂) ; (géométrie B) la fracture F₂ et les diaclases qui l’intersectent, à savoir : R_9a, J₂ et D₆ ; (géométrie C) toutes les fractures identifiées.

Fig. 4 Géométries modélisées par le code 3DEC : (a) géométrie A (fracture F2 seule) ; (b) géométrie B (fractures F2, R9a, J2, D6) ; (c) géométrie C (toutes les fractures – cas de référence). 3DEC model geometries: (a) Geometry A (fracture F2 only); (b) Geometry B (fractures F2, R9a, J2, D6); (c) Geometry C (all the fractures – reference case).

3.1.1 Description du modèle 3DEC

Les simulations ont été réalisées en utilisant la version 4.10 de logiciel 3DEC (Itasca Consulting Group Inc., 1994). La matrice rocheuse et les fractures sont supposées avoir un comportement élastique. La variation de la rigidité de la fracture est supposée être négligeable au cours de l’essai (en raison de la faible amplitude du déplacement mesuré).

3DEC permet de simuler le comportement hydromécanique d’une matrice rocheuse déformable et imperméable découpée par des fractures supposées saturées en eau. L’écoulement dans les fractures est considéré comme laminaire, et il est régi par la loi cubique (Witherspoon et al., 1980). Les calculs hydromécaniques supposent une relation entre les variations d’ouverture hydrauliques de fracture (a − a₀) et la variation de fermeture mécanique (Δu_n), où a₀ est l’ouverture hydraulique à l’état initial avant l’essai, et où la variation de fermeture mécanique (Δu_n) est obtenue en divisant la variation de contrainte effective normale (de Δσ_n′) par la raideur normale de la fracture (kN). La variation de la contrainte normale effective est obtenue en soustrayant la variation de la pression hydraulique à la variation de la contrainte normale totale, selon la loi de Terzaghi (1923).

Les valeurs des données d’entrée du modèle ont été ajustées de sorte à restituer au mieux les mesures de déplacement et de pression :

• –

  propriétés de la matrice rocheuse : densité (ρ) = 2400 kg/m³, module d’Young (E) = 40 GPa, le coefficient de Poisson (ν) = 0,3 ;

• –

  propriétés de la fracture : rigidité normale (k_n) = 5 GPa/m pour les fracture sub-verticale et 50 GPa/m pour les joints de stratification (raideurs tangentielles (k_s) dix fois plus faibles), ouvertures hydraulique (a₀) = 5 × 10⁻⁴ m pour la fracture de la sub-verticale et 10⁻⁴ m pour les joints de stratification.

3.1.1.1 État hydraulique initial

Lorsque la vanne est fermée, un état d’équilibre hydromécanique s’établit dans le versant. Le niveau d’eau remonte et se stabilise environ 7 à 10 m (suivant la distance au versant) au-dessus de la vanne, soit environ 4 m au-dessus du point d’injection en HM₁. Le « petit » modèle, utilisé ici pour interpréter les essais d’injection/pompage, peut donc être considéré comme initialement totalement saturé. Les pressions hydrauliques initiales imposées dans le modèle sont cohérentes avec les valeurs mesurées (45 kPa en HM₁ et 33 kPa en HM₂).

3.1.1.2 État mécanique initial

La répartition des contraintes initiales est plus complexe à estimer car elle dépend notamment de la géométrie du versant et de sa fracturation. Faute de mieux, on a choisi de l’estimer par un calcul préalable en 2D à grande échelle du type de celui présenté dans Guglielmi et al. (2008). Les variations spatiales de l’état de contrainte dans le grand modèle prennent en compte les effets gravitaires et la proximité au versant (nous n’avons pas pris en compte d’éventuels effets liés à la tectonique locale). Le traitement des résultats de ce calcul a permis de déterminer le champ de contrainte initiale en tout point du « petit » modèle comme une fonction de la position de ce point dans l’espace. On considère en particulier que la contrainte totale initiale normale à la fracture instrumentée (F₂) vaut environ 120 kPa en HM₁ et 93 kPa en HM₂. La contrainte effective initiale normale à cette fracture vaudrait donc, selon la relation de Terzaghi, 120 − 45 = 75 kPa en HM₁ et 93 − 33 = 60 kPa en HM₂.

3.1.1.3 Conditions aux limites

On a considéré des conditions aux limites de type « flux nul » en hydraulique et de type « contraintes imposées » en mécanique (les valeurs imposées étant égales aux valeurs initiales résultant du modèle à grande échelle). On a pu établir, que pour la taille choisie du modèle, les conditions aux limites n’avaient pas d’effet significatif sur les résultats.

3.1.1.4 Sollicitations imposées

Elles consistent en une succession de valeurs de pression (essais d’injection) ou de débit (essais de pompage) au point d’intersection du sondage HM₁ avec la fracture F₂. Les simulations hydrauliques et hydromécaniques sont une succession de calculs en régime permanent (un calcul par palier de débit ou de pression). Pour atteindre l’équilibre, 3DEC réalise un certain nombre de cycles hydrauliques. Dans le cas des calculs hydromécaniques, un certain nombre de cycles mécaniques (dits de relaxation), lié à la méthode explicite, est également réalisé afin d’obtenir un équilibre mécanique à chaque cycle hydraulique.

3.1.2 Résultats du modèle de référence

Les variations calculées aux points HM₁ et HM₂ pour le cas de référence (géométrie C) sont donnés dans le Tableau 2. Une comparaison des valeurs calculées avec les mesures montre que le modèle prédit correctement le débit en HM₁ lorsqu’une pression hydraulique est imposée, et réciproquement la pression en HM₁ lorsqu’un débit est imposé (écart inférieur à 12 %). Les déplacements normaux en HM₁ sont également prédits de façon satisfaisant dans le cas des essais d’injection à pression imposée (écart inférieur à 20 %). Sur ce point, les résultats sont moins satisfaisants dans le cas des essais de pompage à débit imposé (écart de 80 %). Cela pourrait être lié à la non prise en compte dans le modèle de la diminution de la raideur normale de la fracture (kn) qui devrait accompagner la diminution de la contrainte normale effective, comme cela a été mise en évidence en laboratoire à partir d’échantillon prélevé sur le site (Lopez et al., 2007 ; Thoraval, 2005). Les variations de pression hydraulique en HM₂ sont correctement prédites dans le cas des essais de pompage à débit imposé (écart inférieur à 13 %). Toutefois, ces pressions sont fortement sous-estimées dans le cas des essais d’injection à pression imposée. Cela pourrait être lié à un effet de la condition à la limite hydraulique (pression supposée constante) qui, en raison de la taille assez réduite du modèle, pourrait réduire artificiellement la variation de pression en HM₂.

Si les résultats obtenus sont acceptables sur certain point, il apparaît donc qu’il faudrait vraisemblablement reconsidérer certaines hypothèses du modèle (taille réduite du modèle, condition aux limites hydrauliques, raideurs des fractures supposée constante) pour obtenir des résultats totalement satisfaisants. Compte tenu des objectifs plutôt méthodologiques de cet article (et de l’importance des temps de calcul), nous avons toutefois estimé que le modèle pouvait être utilisé dans l’état pour une première évaluation des effets des simplifications géométriques (non prise en compte de certaine fracture, homogénéisation). Il reste à vérifier l’impact de ces hypothèses sur les conclusions, notamment dans le cadre de sollicitation hydraulique de plus longue durée (la durée total des expérimentations ne dépasse pas ici quelques heures – voir Sect. 2.3) où la taille du modèle devrait être augmenté.

3.1.3 Effet d’une réduction du nombre de fractures dans le modèle 3DEC

Trois simulations ont été réalisées selon les géométries décrites sur la Figure 4 : géométrie A (fracture F₂ seule), géométrie B (fractures F₂, R_9a, J₂, D₆), géométrie C (toutes les fractures – cas de référence). La Figure 5 montre les iso-valeurs de pressions hydrauliques (a, b et c) et de déplacement normaux (d, e et f) dans le plan de la fracture F₂. On constate logiquement une diminution de la pression hydraulique maximale lorsque la fracture F₂ est connectée à d’autres fractures. Les fractures environnantes ont également un impact sur les déplacements normaux de la fracture F₂ ; on voit clairement, en effet (Figs. 4e et 4f) des discontinuités de déplacement de part et d’autre des traces des fractures complémentaires dans le plan de la fracture F₂.

On a comparé les valeurs de pressions hydrauliques, de déplacements normaux et de débits aux points HM₁ et HM₂ (Fig. 2). On précise ci-après les écarts de prédictions obtenus avec la géométrie « une seule fracture » (géométrie A) relativement à la géométrie « toutes les fractures » (géométrie C) :

• –

  lorsque l’on impose la pression au sein de la fracture F₂ (en HM₁), on calcule les écarts suivants : débit en HM₁ (−15 %) ; pression en HM₂ (+37 %) ; déplacements normaux en HM₁ (−3 %) et en HM₂ (+23 %) ;

• –

  lorsque l’on impose un débit au sein de la fracture F₂ (en HM₁), on calcule les écarts suivants : pression en HM₁ (+16 %) et en HM₂ (+60 %). déplacements normaux en HM₁ (−13 %) et en HM₂ (−40 %).

Les différences entre les géométries B et C sont moins importantes que les différences entre les géométries A et C. En effet, pour le débit, les variations de pression hydraulique et les variations de déplacement normal, l’écart reste inférieur à 3 % en HM₁ et à 15 % en HM₂.

Ces données montrent qu’une simplification excessive de la géométrie du modèle (en ne considérant que la fracture principale instrumentée F₂) peut avoir un impact important sur la réponse hydromécanique et induire des écarts de prédiction importants.

Fig. 5 Iso-valeur de pressions hydrauliques (à gauche) et de déplacements normaux (à droite) dans le plan de la fracture F2 pour les différentes géométries 3DEC : (a) géométrie A (fracture F2 seule) ; (b) géométrie B (fractures F2, R9a, J2, D6) ; (c) géométrie C (toutes les fractures – cas de référence). Contours of hydraulic pressure (left) and normal displacement (right) calculated within the fracture F2 for the different 3DEC model geometries: (a) geometry A (fracture F2 only); (b) geometry B (fractures F2, R9a, J2, D6); (c) geometry C (all the fractures – reference case).

3.2 Développement d’un modèle hydromécanique à double porosité avec FLAC3D

Comme on vient de le voir, les premiers modèles réalisés (à l’aide du code de calcul 3DEC) se limitaient à représenter le comportement hydromécanique des fractures principales du site de Coaraze en négligeant les écoulements à travers les fractures secondaires et la matrice calcaire poreuse, les « blocs » découpés par les fractures principales retenues étant supposés imperméables. Ces modèles assimilaient également le comportement mécanique de ces « blocs » à celui de la matrice rocheuse intacte.

La nouvelle approche vise à associer à chaque « bloc » découpé par les fractures principales retenues, des caractéristiques hydromécaniques équivalentes susceptibles de restituer plus fidèlement le comportement hydromécanique du milieu réel (poreux et fracturé, la fracturation étant liée aux fractures non introduites explicitement) qu’il est supposé représenter (Fig. 6).

La détermination des caractéristiques hydromécaniques équivalentes n’est qu’un préalable aux nouvelles modélisations des expérimentations locales. Celles-ci ont nécessité le choix d’un nouveau code de calcul permettant simultanément la modélisation explicite du comportement hydromécanique des fractures principales (comme le faisait déjà le logiciel 3DEC) et celui des milieux poreux « équivalent », déterminé préalablement.

Fig. 6 Illustration de la méthodologie mise en œuvre pour homogénéiser les propriétés des blocs découpés par les fractures “clés”. Illustration of the upscaling methodology applied to the blocks cut by the “key” fractures.

3.2.1 Description du modèle FLAC3D

Nous avons choisi d’utiliser le logiciel FLAC3D dans sa version 3.10 (Itasca Consulting Group Inc., 2002). Avec ce logiciel, les fractures sont représentées comme des zones tabulaires de roche poreuse de faibles épaisseurs dont les caractéristiques hydromécaniques (déformabilités, perméabilité) sont déterminées analytiquement à partir des caractéristiques des fractures réelles (raideurs, ouverture hydraulique).

Cette démarche permet de transformer les modèles 3DEC précédents en des modèles FLAC3D (Fig. 7) donnant la même réponse hydromécanique à l’injection (dans la mesure où on considère que les zones « interfractures » sont, dans les deux cas, constituées de la matrice calcaire supposée imperméable).

Fig. 7 Construction de la géométrie du modèle FLAC3D à partir du modèle 3DEC : (a) sans homogénéisation ; (b) avec homogénéisation. Construction of the FLAC3D model geometry from the 3DEC model: (a) without upscaling; (b) with upscaling.

3.2.1.1 Détermination des propriétés équivalentes des zones tabulaires représentant les fractures

La procédure a consisté à associer un modèle transversalement isotrope aux zones tabulaires, le plan d’isotropie (1–2) correspondant au plan de la fracture. Les paramètres à spécifier, en plus de l’orientation du plan d’isotropie sont E₁, E₃, G₁₃, ν₁₂ et ν₁₃. Les relations entre ces paramètres et les raideurs normales et de cisaillement (k_n et k_s) sont les suivantes : E₁ = 1.K_s ; E₃ = e.K_n ; G₁₃ = E₁.E₃/(E₁ + E₃) ; ν₁₂ = ν₁₃ = 0 (où e est l’épaisseur de la zone tabulaire).

Par ailleurs, la valeur de la perméabilité intrinsèque de la zone tabulaire) équivalente à celle d’une fracture d’ouverture hydraulique a₀ est : a₀ = k_tabulaire³/12e. On a également estimé la porosité équivalente (n_tabulaire = a/e), ainsi que le module de Biot équivalent (M_tabulaire = K_w/n_tabulaire = e K_w/a, où K_w est la compressibilité de l’eau).

3.2.1.2 Détermination des caractéristiques hydromécaniques équivalentes des blocs fracturés et poreux définis par les fractures principales

Une des techniques d’homogénéisation la plus connue est celle développée par Oda (1986) et améliorée depuis par Ababou (1991). Elle permet de prendre en compte non seulement l’ouverture et l’orientation de chacune des fractures, mais aussi leur extension, ainsi que leurs caractéristiques mécaniques. On détermine alors, pour un volume élémentaire représentatif (VER), le tenseur de conductivité hydraulique équivalent et les tenseurs définissant la relation entre les contraintes et les déformations.

Pour tenir compte des interactions des fractures, d’autres techniques basées sur la simulation numérique par la méthode des éléments finis (Coste et al., 1999 ou par la méthode des éléments distincts. Nous avons mis en œuvre ici celle développée par Thoraval et Renaud (2004). Il s’agit tout d’abord d’identifier la géométrie de chacun des « blocs » à homogénéiser (matrice découpée par des fractures secondaires). On calcule ensuite, par simulations numériques, la réponse de ces « blocs » à différents chargements hydromécaniques, de sorte à déterminer les tenseurs de perméabilité Kij (tenseur d’ordre 2 reliant le vecteur débit et le gradient hydraulique) et de rigidité Tijkl (tenseur d’ordre 4 reliant le tenseur des contraintes au tenseur des déformations) d’un milieu poreux continu « équivalent » au milieu fracturé considéré.

Le Tableau 3 montre les résultats de l’homogénéisation pour les trois géométries décrites précédemment. Logiquement, les propriétés équivalentes varient selon le nombre de fractures explicitement incluses dans les modèles.

3.2.2 Impact de la perméabilité de la roche intacte et de l’homogénéisation sur la qualité des prédictions

On a testé tout d’abord la sensibilité des résultats à la valeur de la perméabilité intrinsèque de la roche intacte. Les résultats montrent que les écarts, par rapport à un modèle supposant une matrice imperméable, restent inférieurs à 1 % aussi longtemps que la perméabilité matricielle reste inférieure à 10⁻¹⁵ m². Si la perméabilité de la masse rocheuse augmente jusqu’à 10⁻¹³ m², les écarts atteignent 15 %.

Le Tableau 4 rassemble les valeurs de pressions hydrauliques et de déplacements normaux calculé en HM₁ et HM₂ dans la fracture instrumentée pour les trois géométries A, B et C considérées. Deux cas de figure sont envisagés selon que les blocs découpés par les fractures principales sont considérés imperméables (comme précédemment avec 3DEC) ou pas. Dans le deuxième cas, les propriétés hydromécaniques « équivalentes » du milieu constitué de la matrice poreuse et des fractures secondaires sont celles données dans le Tableau 3.

On note en particulier que :

• –

  la prise en compte des fractures secondaires par homogénéisation tend à réduire (globalement) les écarts entre la réponse hydromécanique calculée pour les différentes géométries, l’homogénéisation compensant l’absence de ces fractures secondaires. En effet, les écarts calculés entre la géométrie A (une seule fracture principale) et B (quatre fractures principales) passent en moyenne de 19 à 4 % pour le débit en HM₁, de 12 à 5 % pour le déplacement normal en HM₁, de 20 à 10 % pour la pression hydraulique en HM₂. Les écarts entre géométrie A et B restent toutefois importants avec homogénéisation pour ce qui concerne le déplacement normal en HM₂. Dans ce cas, il semble difficile d’estimer correctement le déplacement normal en HM₂ à partir de la géométrie A, même avec homogénéisation des fractures secondaires ;

• –

  la prise en compte, par l’homogénéisation, des fractures secondaires ne semble nécessaire que dans le cas où seule la fracture instrumentée est retenue explicitement (géométrie A). Dans ce cas, la perméabilité équivalente calculée (Tab. 3) est suffisamment importante (environ 2,10⁻⁵ m/s) pour induire une diminution des pressions plus rapide autour du point HM₁ ainsi qu’une augmentation du débit d’injection en ce point. On note sur la Figure 8 une baisse importante des pressions hydrauliques et des déplacements normaux autour le point HM₁ ;

• –

  dès lors qu’un nombre minimal de fractures est retenu explicitement (géométrie B), il n’apparaît plus aussi indispensable de procéder à une homogénéisation. En effet la réponse hydromécanique de la fracture instrumentée en HM₁ et HM₂ est à peu près la même que l’on homogénéise ou non. La comparaison des Figures 8b et 9b confirme que les valeurs de pressions hydrauliques et de déplacements normaux dans la fracture F₂ ne sont pas affectés par le nombre de fractures prises en compte explicitement dès lors que les fractures « exclues » le sont par homogénéisation. Cela s’explique par le fait que, pour ces modalités, la perméabilité équivalente demeure faible (environ 5,10⁻⁷ m/s) devant celle des fractures (environ 10⁻³ m/s pour les diaclases et 10⁻⁵ m/s pour les joints de stratification). Une diminution de sa valeur n’engendre donc pas d’impact sur la réponse hydromécanique de la fracture instrumentée F₂. Les mêmes résultats sont obtenus pour la réponse mécanique de la fracture ; inclure quelques fractures supplémentaires (géométrie B ou lieu de A) réduit le déplacement normal maximum en HM₁ de 17 % sans homogénéisation et de 6 % avec.

Fig. 8 Iso-valeurs des pressions hydrauliques (à gauche) et des déplacements normaux (à droite) calculés au sein de la fracture F2 avec FLAC3D dans le cas de la géométrie A (où seule la fracture F2 est prise en compte de façon explicite) : (a) sans homogénéisation ; (b) avec homogénéisation. Contours of hydraulic pressure (left) and normal displacement (right) calculated within the fracture F2 with FLAC3D in the case of geometry A (where only the fracture F2 is taken into account explicitly): (a) without upscaling; (b) with upscaling.

Fig. 9 Iso-valeurs des pressions hydrauliques (à gauche) et des déplacements normaux (à droite) calculés au sein de la fracture F2 avec FLAC3D en considérant une homogénéisation des fractures non prises en compte de façon explicite : (a) dans le cas de la géométrie B (fractures F2, R9a, J2 et D6) ; (b) dans le cas de la géométrie C (toutes les fractures). Contours of hydraulic pressure (left) and normal displacement (right) calculated within the fracture F2 with FLAC3D considering an upcaling of the fractures not taken into account explicitly: (a) in the case of geometry B (fractures F2, R9a, J2 and D6); (b) in the case of geometry C (all the fractures).

4 Discussion

Les données d’entrée du modèle comprennent des données géologiques, géométriques, mécaniques et hydrauliques. Pour cet article, elles sont spécifiques aux caractéristiques du site de Coaraze. Il est légitime de se demander si les résultats obtenus ici sont extrapolable pour prédire le comportement hydromécanique d’autres sites, si la géologie ou le contexte (pentes rocheuses instables, carrières, stockages souterrains) diffèrent.

Dans le cas du site de Coaraze, les discontinuités majeures étaient les fractures subverticales plus perméables et déformables. Ce n’est pas le cas de tout massif rocheux et cela doit évidemment être considéré pour d’autres applications. Par exemple, dans le cas de fronts rocheux instables ou des carrières à ciel ouvert, les discontinuités majeures pourraient être les plans de glissements potentiels ; dans le cas des stockages, il pourrait s’agir des chemins d’écoulements préférentiels. Il est clair que dans chaque cas, c’est à l’expert d’identifier ces discontinuités majeures en fonction du contexte et des objectifs fixés au modèle prédictif.

À noter que la caractérisation géologique du site de Coaraze a bénéficié d’une nouvelle méthode de relevé de la fracturation par mesures tachéométriques qui présente l’avantage de localiser précisément les fractures dans l’espace. Cette technique apparaît comme complémentaire à la méthode classique par ligne de mesures qui permet uniquement de déterminer les caractéristiques moyennes des différentes familles de fractures (orientation, espacement ou densité). Compte tenu de la sensibilité constatée des prédictions à la prise en compte correcte des fractures « clés », il semble donc judicieux de réserver la nouvelle méthode de mesure à la localisation précise de ces discontinuités.

Les modélisations réalisées sont basées sur des modèles conceptuels représentant le comportement hydraulique, mécanique et hydromécanique des fractures et de la matrice rocheuse poreuse qu’elles découpent. À chacun de ces modèles correspond un certain nombre de paramètres qui peuvent être évalués à l’aide d’essais de laboratoire ou d’expérimentation in situ. Il semble pertinent de se limiter à une caractérisation « classique » pour la matrice rocheuse et les principales familles de fractures. En fonction de la qualité de prédiction recherchée, une caractérisation plus précise et plus complète pourrait être nécessaire pour les discontinuités majeures. À cette fin, des essais in situ du type de ceux réalisés à Coaraze (mesures synchrones des variations de pression et de déformation induite par un essai d’injection/pompage à l’aide d’une sonde amovible) pourraient déterminer avec précision leur perméabilité et leur raideur, ainsi que les termes de couplage éventuels (qui pourraient être important dans le cas des stockages profonds par exemple).

Avant de construire un modèle, il est nécessaire d’identifier clairement l’objectif qu’on lui assigne et le contexte dans lequel on sera amené à l’utiliser. En effet, le degré de complexité que l’on y mettra (ou à l’inverse les simplifications que l’on s’autorisera) dépendent de la possibilité de négliger, de simplifier ou d’homogénéiser tel ou tel phénomène sans que l’on puisse remettre en cause, compte tenu de l’objectif et du contexte de l’étude, la pertinence des résultats obtenus.

Si l’on s’intéresse à la prédiction de la stabilité ou des écoulements à très grande échelle, il n’est pas envisageable de prendre en compte toute la fracturation de manière explicite (pour des problèmes liés à la reconnaissance in situ de cette fracturation et aux limites des outils de modélisation). Il n’est d’ailleurs pas nécessaire de le faire. En effet, cette recherche a montré que le recourt à la technique d’homogénéisation permet de calculer les caractéristiques hydromécaniques d’un milieu continu poreux équivalent, ce qui rend possible une prédiction correcte du comportement hydromécanique global sans introduire explicitement toute la fracturation.

La technique utilisée dans ce programme, basée sur la modélisation de la réponse de la zone à homogénéiser à différents chargements hydromécaniques, peut être remplacée, si on le souhaite, par l’application des formulations analytiques d’Oda (1986), notamment dans le cas où la fracturation que l’on souhaite homogénéiser est connue en moyenne (ou statistiquement).

Le choix de l’outil numérique sera basé sur sa capacité à prendre en compte la géométrie (simplifiée ou homogénéisée), ainsi que les phénomènes physiques retenus après analyse du contexte et des objectifs. Un facteur limitant important est la nécessité, ou non, d’un modèle en trois dimensions. Les codes utilisés dans le cadre de cette recherche couvrent les différents cas de figure. Dans le cas où le milieu rocheux fracturé peut être entièrement homogénéisé, tout code permettant de simuler le comportement hydromécanique des milieux poreux peut convenir. Dans le cas où le milieu rocheux fracturé peut être remplacé par un réseau de fractures découpant une matrice imperméable, tout code permettant de simuler le comportement hydromécanique des fractures peut convenir. En revanche, dans le cas, le plus complexe, où le milieu rocheux fracturé est remplacé par un milieu à double porosité constitué de fractures majeures prises en compte explicitement et d’un milieu poreux équivalent représentant la matrice et l’ensemble des autres fractures, les choix se réduisent. Dans le cadre de cette recherche nous avons adapté le code FLAC3D à ce besoin.

5 Conclusion

À grande échelle, le milieu rocheux fracture est souvent représenté par un milieu continu équivalent car les fracture sont trop nombreuses pour être représenté explicitement dans un modèle informatique. Des techniques appropriées d’homogénéisation doivent être utilisées pour calculer les propriétés hydromécaniques d’un milieu poreux équivalent de sorte d’être en mesure de simuler correctement la réponse hydromécanique du milieu à une sollicitation.

Dans la présente étude, nous avons examiné les simplifications possibles permettant de prédire avec précision son comportement hydromécanique à différentes échelles.

Les principaux résultats de notre étude peuvent être résumés comme suit :

• –

  une approche numérique 3D, couplant 3DEC et FLAC3D, a été développée. Elle permet de simuler le comportement hydromécanique d’un milieu constitué d’un certain nombre de fractures considérées comme majeures, mais également de celui de la roche intacte ou d’un milieu équivalent constitué de la matrice et des fractures d’importance secondaire non introduites explicitement ;

• –

  la réponse hydromécanique in situ d’une fracture à une sollicitation directe (injection/pompage) ne dépend pas que des caractéristiques propres de cette fracture, mais également des caractéristiques de la matrice et des fractures environnantes. La prédiction du comportement de la fracture sollicitée nécessite la prise en compte de l’environnement de la fracture soit explicitement, soit par homogénéisation ;

• –

  la prise en compte explicite de quelques fractures supplémentaires connectées à la fracture sollicitée permet de simuler convenablement les expérimentations locales de quelques heures réalisées dans le site de Coaraze. Il apparaît dans ce cas que la prise en compte des fractures écartées par homogénéisation ne modifie pas significativement les résultats dès lors que le nombre de fractures introduites est suffisant. Cette conclusion pourrait toutefois être mise en défaut dans le cas d’expérimentations de longues durées (plusieurs jours).

Références

Liste des tableaux

Liste des figures

	Fig. 1 Schématisation du site expérimental de Coaraze (Guglielmi, 1999). Description of Coaraze experimental site (Guglielmi, 1999).
Dans le texte

	Fig. 2 Description des expérimentations d’injection et de pompage au sein de la fracture verticale F2. Description of the experimental device used for local injections or pumping within fracture F2.
Dans le texte

	Fig. 3 Sonde de mesure mise en place dans les forages pour mesurer la pression hydraulique et les déplacements normaux dans la faille F2. Experimental device inserted into the boreholes to measure hydraulic pressure and normal displacements within the fault F2.
Dans le texte

	Fig. 4 Géométries modélisées par le code 3DEC : (a) géométrie A (fracture F2 seule) ; (b) géométrie B (fractures F2, R9a, J2, D6) ; (c) géométrie C (toutes les fractures – cas de référence). 3DEC model geometries: (a) Geometry A (fracture F2 only); (b) Geometry B (fractures F2, R9a, J2, D6); (c) Geometry C (all the fractures – reference case).
Dans le texte

Fig. 5 Iso-valeur de pressions hydrauliques (à gauche) et de déplacements normaux (à droite) dans le plan de la fracture F2 pour les différentes géométries 3DEC : (a) géométrie A (fracture F2 seule) ; (b) géométrie B (fractures F2, R9a, J2, D6) ; (c) géométrie C (toutes les fractures – cas de référence). Contours of hydraulic pressure (left) and normal displacement (right) calculated within the fracture F2 for the different 3DEC model geometries: (a) geometry A (fracture F2 only); (b) geometry B (fractures F2, R9a, J2, D6); (c) geometry C (all the fractures – reference case).

Dans le texte

	Fig. 6 Illustration de la méthodologie mise en œuvre pour homogénéiser les propriétés des blocs découpés par les fractures “clés”. Illustration of the upscaling methodology applied to the blocks cut by the “key” fractures.
Dans le texte

	Fig. 7 Construction de la géométrie du modèle FLAC3D à partir du modèle 3DEC : (a) sans homogénéisation ; (b) avec homogénéisation. Construction of the FLAC3D model geometry from the 3DEC model: (a) without upscaling; (b) with upscaling.
Dans le texte

Fig. 8 Iso-valeurs des pressions hydrauliques (à gauche) et des déplacements normaux (à droite) calculés au sein de la fracture F2 avec FLAC3D dans le cas de la géométrie A (où seule la fracture F2 est prise en compte de façon explicite) : (a) sans homogénéisation ; (b) avec homogénéisation. Contours of hydraulic pressure (left) and normal displacement (right) calculated within the fracture F2 with FLAC3D in the case of geometry A (where only the fracture F2 is taken into account explicitly): (a) without upscaling; (b) with upscaling.

Dans le texte

Fig. 9 Iso-valeurs des pressions hydrauliques (à gauche) et des déplacements normaux (à droite) calculés au sein de la fracture F2 avec FLAC3D en considérant une homogénéisation des fractures non prises en compte de façon explicite : (a) dans le cas de la géométrie B (fractures F2, R9a, J2 et D6) ; (b) dans le cas de la géométrie C (toutes les fractures). Contours of hydraulic pressure (left) and normal displacement (right) calculated within the fracture F2 with FLAC3D considering an upcaling of the fractures not taken into account explicitly: (a) in the case of geometry B (fractures F2, R9a, J2 and D6); (b) in the case of geometry C (all the fractures).

Dans le texte