<figure><img src="image_1.png" /></figure>

HAL
open science

# Modelling the creep behaviour and induced failure of clay rock from microscale viscosity to large-scale time-dependant gallery convergences using a multiscale numerical approach

Yufeng Sun, Benoît Pardoen, Henry K.K. Wong

► To cite this version:

Yufeng Sun, Benoît Pardoen, Henry K.K. Wong. Modelling the creep behaviour and induced failure of clay rock from microscale viscosity to large-scale time-dependant gallery convergences using a multiscale numerical approach. Computers and Geotechnics, 2023, 162, pp.105691. 10.1016/j.compgeo.2023.105691. hal-04236855

HAL Id: hal-04236855

https://hal.science/hal-04236855v1

Submitted on 1 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Modelling the creep behaviour and induced failure of clay rock from microscale viscosity to large-scale time-dependent gallery convergences using a multiscale numerical approach

Yufeng SUN*, Benoit PARDOEN, Henry KK WONG

Université de Lyon, ENTPE, LTDS, Vaulx-En-Velin, France

* Corresponding author: yufeng.sun@entpe.fr

7

## 8 Abstract

The large-scale creep behaviour of Callovo-Oxfordian (COx) clay rock is modelled from small-scale viscous mechanisms of the rock using a multiscale numerical approach in the context of radioactive waste repositories. At the mesoscale, the saturated non-homogeneous rock is represented in digital 2D Representative Elementary Areas (REAs) as a cracked composite material consisting of rigid elastic mineral inclusions (quartz, calcite, and pyrite) embedded in a clay matrix. The modelling of these materials is considered within double-scale finite element framework, in which the homogenised responses of the REA to the enforced global kinematics serve as a numerical constitutive relation at the macroscale. To reproduce the rock creeping, two viscous mechanisms have been introduced to consider the creep of the clay matrix: either the viscoplasticity of the clay aggregates or the viscoelasticity of their contacts, or both. Firstly, laboratory biaxial creep tests of clay rock samples are simulated. A three-stage creep stage is reproduced and the creep failure process is discussed. Then, large-scale underground engineering structures are modelled to reproduce its time-dependent behaviour. It is found that the developed multiscale model is able to provide some valuable insights into the large-scale creep behaviour of clay rocks through the morphological and material small-scale characterization of REA.

**Key words:** Nuclear waste disposal; COx clay rock; Double scale FE²; Homogenised response; Creep
24

# 1. Introduction

25

26

27

28

29

30

31

32

33

34

35

36

37

The French National Radioactive Waste Management Agency (Andra) has built the Meuse/Haute-Marne Underground Research Laboratory (MHM URL) at Bure (East of Paris sedimentary basin) since 2000 aiming to demonstrate the feasibility of an industrial radioactive waste repository in deep geological formation. The target formation is the Callovo-Oxfordian (COx) clay rock, with depths between 420 and 550 m below the ground surface, with properties that are very favourable for radioactive waste disposal, such as their low hydraulic permeability (10<sup>-22</sup> to 10<sup>-20</sup> m<sup>2</sup>), good resealing capacity, and self-sealing properties. From the mineralogy study and microscopic observations, the COx clay rock is mainly composed of four components: tectosilicates (–Si<sub>x</sub>O<sub>y</sub>, mainly quartz, 10–40%), carbonates (–CO<sub>3</sub>, mainly calcite, 15–80%), heavy minerals (FeS<sub>2</sub>, pyrite, in a low proportion of 0–3%), and clay minerals (20–60%) (Armand et al., 2017). The URL is built at an average depth of 490 m corresponding to the median depth of the Callovo-Oxfordian geological formation. At this level, in the clay-rich lithostratigraphic unit of the COx clay rock, the proportions of quartz, calcite, pyrite and clay minerals are respectively 20%, 20–25%, 3% and 50–55% (Armand et al., 2014).

38

39

40

41

Various studies undertaken by Andra include the development of theoretical and numerical models to characterise the complex COx clay rock behaviour in the short- and long-terms. This requires the determination of constitutive models and parameters based on experimental results, then the numerical simulation of underground galleries including construction and exploitation stages (Seyedi et al., 2017). The long-term creep deformation and viscous behaviour of the COx clay rock is one of the key issues being investigated in underground experiments, to ensure the safety conditions required for long-term repository of radioactive wastes. In many cases, the creep behaviour of the clay rock is described by using phenomenological viscoplastic or viscoelastic models (e.g. Pardoen and Collin, 2017; Sun et al., 2021);

42

43

44

45
however, these models cannot give a clear explanation on the influence of microstructural properties on the material behaviour. The complex multi-field coupling phenomenon of underground structures requires realistic modelling and to give more reasonable insight of structural behaviour from small-scale characteristics (Fu et al., 2020; Yuan et al., 2021). From this perspective, a material constitutive model should detail the material small-scale structure (mesoscopic, microscopic, etc.) relevant to reproduce the macroscale material behaviour. Theoretically and numerically, this involves separate descriptions of each material constituent and their interactions at microscale, the constituent spatial arrangement and heterogeneity at mesoscale, then scale transition towards material behaviour at macroscale. Multiscale numerical analysis such as FEMxFEM or FE² (van den Eijnden et al., 2016), which is suitable for quasi-brittle materials as the considered clay rock, and FEMxDEM (Desrues et al., 2019), which is suitable for granular materials, can effectively model the multi-scale behaviours of geomaterials. These specific methods allow to define a material Representative Element Area (REA) containing the material mesostructure information at each material point (Gaussian integration point) of the macroscale FEM computation. The homogenised mesoscale responses of the REA, obtained from microscale FEM computation, represent the macroscale material behaviour at each material point.

Considering the short-term excavation of underground galleries in quasi-brittle clay rock, the induced rock damage and fracturing processes are important to reproduce (Pardoen et al. 2015b). They induce the development of an Excavation Damaged Zone (EDZ) around underground galleries and nuclear waste cells. Moreover, in presence of underground water, the hydromechanical (HM) coupled behaviour of rock is also important to reproduce, especially in damaged and fractured rock media (Armand et al. 2013; Mánica et al. 2021; Xu et al., 2018; Yang et al., 2023; Zhan et al. 2021). Consequently, microscale models have been developed to reproduce damage and solid-fluid interactions at the mineral grain scale,
as for instance the FE model proposed by Frey et al. (2013). Such model allows to capture the degradation of geological materials, the physical processes of water flow, and the hydromechanical coupling (e.g. hydraulic permeability evolution in microcracks). Further developments and studies have allowed to extend this model and to confirm the capabilities of the double-scale FE² modelling in capturing the characteristic behaviour of a damaged rock at microscale as well as of an excavation fractured zone at macroscale, through qualitative analysis (van den Eijnden et al., 2016, 2017; Mourlas et al., 2023; Pardoen et al., 2020). In this model, the delayed behaviour of the material is induced by the pore water advection process in the crack (van den Eijnden et al., 2016). If the rock porous space is dry there is no longer delayed behaviour, but only elastic strain. Recently, Bertrand et al. (2020) developed a multiscale model for the modelling of multiphase flow in fractured reservoirs, particularly in coal seams during coalbed methane recovery (or carbon dioxide storage). By comparing the consolidation characteristic times predicted by the poro-elastic theory with experimental results on clay rock, Rahal et al. (2017) found that the delayed clay rock behaviour cannot be explained by pore pressure dissipation alone, because it includes creep phenomena due to a high clay content. Therefore, they conclude that damage, poro-mechanics and viscoplasticity have to be coupled if complex behaviour of clay rock around the gallery is to be explained.

A short description of the COx clay rock behaviour at different scales is given in Section 2. The double scale FE² model is introduced in Section 3, with macroscale field equations, mesoscale formulation of the periodic boundary value problem (BVP), and microscale constitutive models. Then, the results of double-scale numerical simulation at laboratory scale are detailed in Section 4 where the different creep phases, including the tertiary creep process, are discussed. Lastly, the results of double-scale numerical simulation at underground gallery scale are discussed in Section 5.
90

## 2. Creep and fracturing behaviour of COx clay rock at different scales

91

92

93

94

95

96

97

98

100

101

At mesoscale, clay rocks are heterogeneous and composed of several mineral inclusion types embedded in a clay matrix (Cosenza et al., 2015a, 2015b; Robinet, 2008; Robinet et al., 2012). 2D and 3D spatial distributions of the mineral groups of the COx clay rock obtained from the analyses of X-ray Computed micro-Tomography (micro-CT) and Scanning Electron Microscopy (Backscattered Electrons, SEM BSE) images are shown in **Fig. 1**. These analyses lead to the determination of mineral groups, spatial distributions, mineral area fractions, and mineral morphology (size, shape, and orientation of mineral inclusions). Therefore, at the mesoscale, the COx clay rock will be considered as a solid four-phase composite material in which three different types of mineral inclusions (quartz, calcite, and pyrite) are embedded in a clay matrix. From the work of Cosenza et al. (2015a, 2015b), the average contents of these four minerals in the COx clay-rich lithostratigraphic unit are 18%, 30%, 2%, and 50% based on detailed experimental investigations.

<figure><img src="image_2.png" /></figure>

102

103

104

105

106

**Fig. 1.** Spatial distribution of mineral groups of the COx clay rock: (a) 2D SEM image on a section perpendicular to the bedding planes; (b) 3D micro-CT subvolume (EST26095) (Robinet et al., 2012).
As an approximation, three relevant scales of the clay rock behaviour are considered and introduced in the background of this paper:

107

* Macroscale: laboratory sample scale, large gallery scale, EDZ.

108

109

Once the micro-damage threshold is reached (at the microscale of mineral contact), microcracks and mesocracks initiate (in and through the mesostructure), then grow, accumulate, and propagate
within the material at macroscale. If distributed microcracks start to coalesce, it further leads to the development of strain localisation (e.g. shear bands for shear strain localisation). To represent the macro-fracturing process, macro shear strain localisation is considered as a precursor to macro-fractures. Concerning the creep behaviour at macroscale, it is considered as coming from clay matrix viscosity at smaller scales and can be influenced by mechanical interactions between a set of REAs.

* Mesoscale: scale of the arrangement of mineral inclusions and of the connected clay matrix which are considered in a REA.

Under loading, the decohesion appears in a localised manner for a small proportion of mineral contacts (Pardoen et al., 2020). If a coalescence of fully damaged mineral contacts occurs, a microcracking path develop between the mineral inclusions and goes through the clay matrix (intergranular cracks). The microcracking propagates then through the whole mesostructure (REA) to form a meso-crack. The overall material behaviour is the homogenised behaviour over the REA.

* Microscale: scale of the mineral inclusions, clay aggregates, fracture porosity.

The micro characteristics include the inclusion morphology and orientations. The local behaviour is considered in each solid phase (elastic mineral inclusion and viscosity in clay aggregates) and at their contact (damageable cohesive interface, fluid flow and viscous sliding). Stress loading can engender micro-damage at the scale of the mineral contacts, by decohesion. Concerning the creep behaviour, the viscosity in clay aggregates and their contacts leads to the viscosity of the clay matrix, which is composed of the assembly of the clay aggregates and of the aggregate contacts.
Examples and evidences of the creep behaviour of the COx clay rock at different scales are depicted in **Fig. 2** (after Armand et al., 2013, 2014; Pardoen and Collin, 2017; Sun et al., 2023a). Time evolution of the convergence of large-scale underground galleries at the MHM URL has been observed for a long period (Armand et al., 2013). For example, the *in situ* monitoring of gallery wall convergence (**Fig. 2 (a)**) for the Galerie Expérimentale Deux (GED) has been carried out by Andra for more than 10 years. No evidence of convergence stabilisation has been observed, but the deformation rate gradually slows down with time. A model benchmark exercise has been launched by Andra since 2012 to provide an overall view of the developed models on the *in situ* observations (Seyedi et al., 2017). Several hydromechanical coupling models (van den Eijnden et al., 2017; Jung et al., 2022; Mánica et al., 2021; Pardoen et al., 2015a, 2015b; Pardoen and Collin, 2017; Souley et al., 2017) have been proposed and aim to capture the long-term hydromechanical characteristics of *in situ* galleries.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Creep behaviour of COx clay rock at different scales.</figcaption></figure>
Then, at the scale of laboratory test, irreversible creep strain evolution with time have been measured on pluricentimetric laboratory COx clay rock samples subjected to constant stress (Fig. 2 (b)). The laboratory results indicate that the creep process of clay rock under constant deviatoric stress can be composed in three stages (Liu et al., 2018; Shahbodagh et al., 2020; Weng et al., 2023): a primary creep in the short term with a decrease of the creep strain rate, a secondary creep with a constant creep strain rate over time (i.e. steady state creep), and a possible tertiary creep in the long term with an increase of the creep strain rate towards creep failure (Fig. 3). The creep characteristics can be observed in different types of rock, and their creep curves are basically similar despite their different strengths (Chen et al., 2018). The increase of the creep rate appears to be related to the onset and development of damage (Liu et al., 2018). This can further lead to a creep-induced failure of the material. Up to now, few results have been obtained so far on tertiary creep for COx clay rocks (Liu et al., 2018; Zhang et al., 2019). Reproducing such complete creep behaviour of clay rock using a multiscale approach is one of the objectives in this paper. In some studies, the existence of a stress threshold from which viscoplastic strains start to develop has been highlighted by multi-step creep tests (Fabre and Pellet, 2006) and confirmed by one-step creep test (Liu et al., 2015). Moreover, the laboratory tests on COx clay rock conducted by Armand et al. (2017) show that the deviatoric loading level and confining pressure both have an important effect on the amplitude of the creep strain.
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Schematic diagram of creep stages and creep failure induced by the material creep with time (adapted from Frenelus et al. (2022)).</figcaption></figure>

At the meso- and micro- scales, the deformation in clay rocks is dominated by damage, decohesion, and cracking mechanisms (Desbois et al., 2017). With a combination of SEM and broad ion beam (BIB) carried out by Desbois et al. (2017), microstructural observations on COx clay rock at microscale show that the potential decohesion mechanisms around mineral inclusions and cracking within the clay matrix develop in the clay rock. The microcracking includes inter-granular microfractures propagating in the clay matrix and between mineral inclusions and, as well as intragranular and trans-granular microfractures propagating in non-clay minerals. REAs have to describe the grains but also the failure and damage modes at small scale. To do so, the interfaces between two different mineral phases (inclusion-clay) and within the clay matrix (clay–clay aggregate contacts) are considered as potential microcracks, as shown in **Fig. 2 (c)**. However, several micromechanisms during deformation process of the COx clay rock such as grains break and pore collapse are not considered in the model. In order to model the creep behaviour of the clay rock from a small-scale mechanics, the viscosity is introduced by the development of two viscous mechanisms, the “diffuse” viscoplastic flow of the clay aggregates or the viscous sliding between large rigid clay aggregates (at their contacts) (Sun et al., 2023b).
179

### 3. Double scale FE² model

180

181

182

183

184

185

186

187

The used double-scale modelling framework is the so-called finite element squared (FE²) method which allows scale transition of the material behaviour by computational homogenisation (Fig. 4), aiming at improving the modelling of hydromechanical behaviour of complex heterogeneous rocks using a multiscale modelling approach. This approach is of particular interest for both underground and computational engineering, aiming to bridge the gap between material and (geo)structural scale modelling. At the same time, several limitations exist in current modelling approach, for example a number of mesostructural and physical parameters that need to be calibrated, time-consuming calculation, and the necessity of a more realistic 3D mesostructure.

188

189

190

191

192

193

194

195

196

197

198

199

As mentioned previously, within the framework of computational homogenisation, a mesoscale BVP is carried out on a REA of which the homogenised responses are considered as numerical constitutive relations at macroscale. The method allows a dialog between the two computations: from macro to meso scale (downscaling), several macroscale variables (e.g. the strain tensor, fluid pressure, fluid pressure gradient) are transferred to the mesoscale BVP as the REA boundary conditions; from meso to macro scale (upscaling), averaging the meso-BVP solution allows to compute the macroscale dual quantities (e.g. the stress tensor, fluid mass flux, fluid mass content rate) and tangent stiffness matrix. Furthermore, the entire material mesostructure is represented and modelled in a REA. Different components inside the REA are described by their individual constitutive models. The framework of computational homogenisation has been successfully extended for multiphysics coupling problems such as hydromechanical coupling (van den Eijnden et al., 2016; Frey et al., 2013; Marinelli et al., 2016), partially saturated media (Bertrand et al., 2020), and thermal coupling (Ozdemir et al., 2008; Zalamea et al., 2021).
200

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Schematic diagram of FE² model (adapted from van den Eijnden et al., 2016).</figcaption></figure>

201

202

203

204

205

206

207

208

209

210

211

212

213

214

As for the mesostructural (REA) modelling, an algorithm based on Voronoï tessellation from the work of Pardoen et al. (2020) and van den Eijnden et al. (2016, 2017) is adopted to generate periodic 2D REAs with random mesostructures. The representativeness of the mesostructure is evaluated by comparison of numerical and experimental geometrical characteristics of the minerals. The characteristic length of the mesostructure is $L_{REA} \approx 100$ µm defined experimentally (Cosenza et al., 2015a, 2015b). For each type of mineral, Pardoen et al. (2020) give a good agreement between the numerical and experimental results: the minimal-mean-maximal sizes are: 6–72–230 µm² for quartz, 6–43–115 µm² for carbonates, and 6–35–84 µm² for pyrite. The morphology of the mineral inclusions is also evaluated by ellipse fitting procedure with average measurements, as summarised in **Table 1** by Pardoen et al. (2020). In order to obtain mesostructures with grains with a certain preferential orientation, van den Eijnden et al. (2017) have introduced a grain shape tensor for stretching and rotating the mesostructure with respect to the periodic frame.

215

216

<table><caption>**Table 1** Mineralogical and microstructural average characteristics of COx clay rock in vertical sections perpendicular to the bedding planes (EST26095).</caption><thead><tr><td>Minerals</td><td>Area fraction (%)</td><td>Area (µm²)</td><td>Elongation (-)</td><td>Orientation (°)</td></tr></thead></table>
<table><tr><td>Tectosilicates (quartz)</td><td>18</td><td>100</td><td>0.64</td><td>0 (variable)</td></tr><tr><td>Carbonates (calcite)</td><td>30</td><td>45</td><td>0.55</td><td>0</td></tr><tr><td>Heavy minerals (pyrite)</td><td>2</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Clay matrix</td><td>50</td><td>-</td><td>-</td><td>-</td></tr></table>

## 3.1. Macroscale formulation

Tensile stresses and strains are considered positive here below, following material mechanic sign convention. The macroscale behaviour of the considered COx clay rock is quasi-brittle under in situ compression stress range and its failure is related to the development of macroscale fractures (Pardoen et al., 2015b). In a double-scale FE² approach, the latter results from the damage and cracking at meso and micro scales (van den Eijnden et al., 2016, 2017; Mourlas et al., 2023). Such macroscale strain softening behaviour is due, in FE continuous modelling, to the appearance of strain localisation zones (Pardoen et al., 2015a, 2015b; Pardoen and Collin, 2017). Consequently, the adopted double-scale approach is combined with a regularisation method, i.e. the local second gradient model (Collin et al., 2006), to avoid macroscale mesh dependency in FE methods involving strain localisation (Fig. 4). To implement this theory into finite element code, two governing equations are introduced through a field of Lagrange multipliers $\lambda_{ij}$:

$$ \int_{\Omega^t} \left( \sigma_{ij}^t \frac{\partial u_i^*}{\partial x_j^t} + \Sigma_{ijk}^t \frac{\partial v_{ij}^*}{\partial x_k^t} - \lambda_{ij}^t \left( \frac{\partial u_i^*}{\partial x_j^t} - v_{ij}^* \right) \right) d\Omega^t = \int_{\Gamma^t_\sigma} (\bar{t}_i^t u_i^* + \bar{T}_i^t v_{ij}^* n_j^t) d\Gamma^t \quad (1) $$

$$ \int_{\Omega^t} \lambda_{ij}^* \left( \frac{\partial u_i^t}{\partial x_j^t} - v_{ij}^t \right) d\Omega^t = 0 \quad (2) $$

Eq. (1) is the weak form of the momentum balance equation between internal and external virtual works at current material configuration (unit volume) $\Omega^t$, in which the quantities with superscript * denote virtual quantities; $\sigma_{ij}^t$ are the components of the Cauchy stress tensor, $x_j^t$ are the current coordinates, $u_i^*$
232 is the macroscale displacement field, $v_{ij}^*$ is the microkinematic gradient field. Furthermore, $\Sigma_{ijk}^t$ are the components of the double (second order) stress tensor. In the right hand side of Eq. (1); $\bar{t}_i^t$ is the classical external traction per unit area; $\bar{T}_i^t$ is an additional external double force per unit area, both applied on the part $\Gamma_o^t$ of the boundary $\Gamma^t$ of $\Omega^t$, and $n_j^t$ is the normal unit vector to the boundary.

233
234
235
236 For the fluid part, a fluid mass balance equation is needed for the virtual fluid pressure $p^*$ and can be written in the weak form as:

$$
\int_{\Omega^t} \left( \dot{M}^t p^* - m_i^t \frac{\partial p^*}{\partial x_i^t} \right) d\Omega^t = \int_{\Omega^t} Q^t p^* d\Omega^t - \int_{\Gamma^t} \bar{q}^t p^* d\Gamma \quad (3)
$$

238 where $M^t$ is the fluid mass with $\dot{M}^t$ its time derivative, $m_i^t$ is the fluid mass flux, $Q^t$ is a sink term and
239 $\bar{q}^t$ is the boundary input flux per unit area which can be calculated by the product of $m_i^t$ and the boundary
240 surface normal vector $n_i^t$.

241 **3.2. Mesoscale formulation**

242
243
244
245
246
247
248
249
250
The COx clay rock is assumed as a cracked heterogeneous medium at mesoscale, in which rigid mineral
grains (mainly quartz, calcite, and pyrite) are embedded within a soft clay matrix. To model failure and
damage modes at small scale, the contact interfaces between solid mineral grains are assumed as
microcracks and correspond to two types of failure: microcracks between two different mineral phases
(inclusion-clay) and within the clay matrix (clay–clay). These interfaces form a network of pore channels
through which fluid can penetrate and flow (Fig. 5). The creep deformations of host rock are reproduced by
introducing viscoplastisity into clay aggregates or viscoelasticity at clay aggregate contacts. The realistic
periodic mesostructure of the COx clay rock suitable for FE² numerical approach has been introduced by
van den Eijnden et al. (2017) and Pardoen et al. (2020), with 2D random Voronoi tessellation, in which
several geometric and physical properties can be reproduced parametrically, for example, granulometry, elongation, angularity, etc.

<figure><img src="image_6.png" /><figcaption>Fig. 5. Schematic diagram of fluid network at microscale in present FE² model.</figcaption></figure>

### 3.2.1. Mesostructure Generation

### 3.2.2. Separation of scales

In order to meet the requirement of macroscale homogenisation, the small strain assumption is assumed to overcome the definition problems of stress and strain states in and around the interfaces at the microscale. The macroscale deformation gradient tensor $\underline{F}$ is decomposed into a rotational component $\underline{R}$ and a symmetric strain component $\underline{U}$ (van den Eijnden et al., 2016):

$$F_{ij} = \frac{\partial x_i}{\partial x_j} = R_{ik}U_{kj} \qquad (4)$$

The small strain tensor enforced on the REA is then defined as $\underline{\epsilon} = \underline{U} - \underline{\delta}$, in which $\underline{\delta}$ is the Kronecker delta tensor. Both kinematics fields at meso and macro scales exist and are different in double-scale FE² framework.
Given a material point $\hat{P}$ of position $\vec{x}$ with a displacement at macroscale $u_i(\hat{x})$ and fluid pressure $p(\hat{x})$, the separation of scales for displacement and fluid pressure fields are formulated as (van den Eijnden et al., 2016):

$$ \frac{\partial u_i(\hat{x})}{\partial x_j}(x_j - \hat{x}_j) + u_i^f(x) \ll u_i(\hat{x}) \qquad (5) $$

$$ \frac{\partial p(\hat{x})}{\partial x_j}(x_j - \hat{x}_j) + p^f(x) \ll p(\hat{x}) \qquad (6) $$

where $\vec{x}$ is the position of a material point close to $\hat{P}$, $u_i^f(x)$ and $p^f(x)$ represent the mesomechanical fluctuation field of displacement and fluid pressure respectively. Eq. (5) means that the length scale of the REA should be much smaller than the macroscale problem. Simply written, $x_i - \hat{x}_i \ll u_i(\hat{x})$.

### 3.2.3. REA with periodic boundary conditions

The macroscale kinematics enters the mesoscale BVP through boundary conditions enforced on the REA external boundaries. For the problems with elliptic equations underlying the BVP of REA, the periodic boundary condition is adopted as it is the most efficient way to enforce global kinematics on the REA (van der Sluis et al., 2000; Terada et al., 2000). The external domain boundary of the REA is the place where periodic boundary conditions are imposed and is subdivided into two parts: the lead part $\Gamma^L$ and the follow part $\Gamma^F$. The kinematics of any point on the follow boundary depend on the kinematics of its corresponding homologous point on the lead boundary. The distance between these two points is defined by the period vector $y_i = x_i^F - x_i^L$.

The mechanical part of the periodic boundary condition of the REA is defined in terms of the displacement relationship between homologous points as:

$$ \vec{u}^F(x^F) = \vec{u}^L(x^L) + \nabla \vec{u} \cdot \vec{y} \qquad (7) $$

In the same way, the periodic boundary condition for hydraulic part writes:
$$p^F = p^L + \nabla p \cdot \vec{y} \qquad (8)$$

Moreover, the boundary traction $\vec{t}$ and boundary fluid mass flux $\bar{q}$ need to satisfy the antiperiodic condition:

$$\vec{t}^F + \vec{t}^L = 0 \qquad (9)$$

$$\bar{q}^F + \bar{q}^L = 0 \qquad (10)$$

### 3.2.4. Homogenised mesoscopic responses

After solving mesoscale BVP as shown in Fig. 4, the homogenised responses including stress tensor $\sigma_{ij}$, fluid mass flux $m_i$, and fluid mass rate $\dot{M}$ are obtained using homogenisation method based on Hill's condition (Hill, 1965; Mandel, 1972). This allows to write the homogenised macroscale responses as:

$$\sigma_{ij} = \frac{1}{2\Omega} \int_{\Gamma^F} (t_i y_j + t_j y_i) d\Gamma \qquad (11)$$

$$m_i = \frac{1}{\Omega} \int_{\Gamma^F} \bar{q}^F y_i d\Gamma \qquad (12)$$

$$\dot{M} = \frac{M^t - M^{t-\Delta t}}{\Delta t} \qquad (13)$$

where $\Delta t$ is current incremental time step. In fact, the fluid mass rate $\dot{M}$ can be obtained from Eq. (13) directly instead of using homogenisation process.

## 3.3. Microscale model

The COx clay rock is assumed to be composed of an assembly of deformable grains connected by cohesive interfaces at the mesoscale (as shown in Fig. 5). Several constitutive models at the microscale are therefore introduced in this section. Specifically, Section 3.3.1 and Section 3.3.2 respectively show the interface
models with and without viscosity, Section 3.3.4 describes the channel flow model of interfaces, and
Section 3.3.3 describes the constitutive relations of solid grains.

### 3.3.1. Interface cohesive model

The contact interfaces between solid grains are modelled as damaged and cohesive, and the interface cohesive models in the normal and tangential directions to the grain boundary are respectively shown in Fig. 6. The cohesive forces at two directions are formulated as (van den Eijnden et al., 2016):

$$c_n^t = \begin{cases} c_n^{\max}(1 - D_n^t) \frac{1}{D_n^t} \frac{\Delta u_n^t}{\delta_n^c}, & \text{if } \Delta u_n^t \ge 0 \\ c_n^{\max}(1 - D_n^t) \frac{1}{D_n^t} \frac{\Delta u_n^t}{\delta_n^c} - \kappa \Delta u_n^{t2}, & \text{if } \Delta u_n^t < 0 \end{cases} \quad (14)$$

$$c_t^t = c_t^{\max}(1 - D_t^t) \frac{1}{D_t^t} \frac{\Delta u_t^t}{\delta_t^c} \quad (15)$$

where subscripts ‘n’ and ‘t’ respectively indicate normal and tangential directions; $c_{n/t}^t$ are the interface cohesive forces at time $t$; $c_{n/t}^{\max}$ are the maximum cohesions; $\Delta u_{n/t}^t$ are the interface relative displacements (relative to the unit REA size) at time $t$; $\delta_{n/t}^c$ are the critical relative displacements for complete decohesion; $\kappa \gg c_n^{\max}$ is a penalisation term to avoid grain inter-penetration. Considering loading-unloading conditions, $D_{n/t}^t$ are the interface damage state parameters which depends on time history ($\tau = 0 \dots t$):

$$D_{t/n}^t = \min \left( \max \left( D_{t/n}^0, \frac{1}{\delta_{t/n}^c} \max \left( |\Delta u_{t/n}^\tau| \right) \right), 1 \right) \quad (16)$$

where $D_{t/n}^0$ is the critical point at which the interface damage begins. Therefore, if $D_{t/n}^t = D_{t/n}^0$, no damage occurs; if $D_{t/n}^0 < D_{t/n}^t < 1$, the interface is in the softening stage; and if $D_{t/n}^t \ge 1$, the interface is in the complete decohesion state.
<figure><img src="image_7.png" /><figcaption>Fig. 6. Damageable cohesive interface model in (a) normal direction; and (b) tangential direction.</figcaption></figure>

### 3.3.2. Viscoelastic clay aggregates contacts

When the viscosity is considered in clay aggregate contacts, the viscoelastic model is used combining with the one described in Section 3.3.1. The schematic representation of viscoelastic model as shown in Fig. 7 is commonly denominated as the Standard Linear Model or alternatively Zener model (Zener, 1948). In this model, the total relative displacement $\Delta u_i^{ve}$ ($i = t, n$) can be written as the sum of its elastic $\Delta u_i$ and viscous parts $\Delta u_i^v$:

$$ \Delta u_i^{ve} = \Delta u_i + \Delta u_i^v \qquad (17) $$

<figure><img src="image_8.png" /><figcaption>Fig. 7. Schematic representation of viscoelastic model of clay aggregate contacts.</figcaption></figure>

Similarly, the total cohesion force $c_i^{ve}$ writes (Sun et al., 2023b):

$$ c_i^{ve} = c_i + c_i^v \quad ; \quad c_i^v = \beta c_i = \beta \frac{c_i^{max}(1-D_i^t) \Delta u_i^{ve} - \Delta u_i^v}{D_i^t \delta_i^c} \qquad (18) $$
where $\beta$ is a material parameter describing the relation between the elastic and the viscous cohesion forces. The evolution law for the internal variable $\Delta\dot{u}_i^v$ is formulated as:

$$ \Delta\dot{u}_i^v = \frac{\delta_i^c c_i^v}{\mu_i} \qquad (19) $$

where $\mu_i$ is the viscosity of interfaces in normal ($i=n$) or tangential ($i=t$) direction. The evolution of viscoelastic opening $u_{ri}^v$ can be expressed in the integration form as:

$$ \Delta u_i^v = \Delta u_i^{ve} - \Delta u_i = \Delta u_i^{ve} - \int_{-\infty}^{t} e^{\frac{t-\tau}{\lambda_i}} \Delta \dot{u}_i^{ve}(\tau) d\tau \qquad (20) $$

where $\lambda_i = \frac{\mu_i}{\left(\frac{c_i^{max}(1-D_i^0)}{D_i^0}\right)}$ is the relaxation time. Assuming in the n-th time step, the times at the start and the end of the increment are respectively denoted as $t_n$ and $t_{n+1}$. After some mathematical calculations, we can get $c_i^v$ at $t_{n+1}$ as follows:

$$ c_i^v = \frac{c_i^{max}(1-D_i)}{D_i \delta_i^c} (\Delta u_i)_{n+1} = \beta \frac{c_i^{max}(1-D_i)}{D_i \delta_i^c} \left( e^{-\frac{\Delta t}{\lambda_i}} (\Delta u_i)_n + \frac{1-e^{-\frac{\Delta t}{\lambda_i}}}{\frac{\Delta t}{\lambda_i}} \delta \Delta u_i^{ve} \right) \qquad (21) $$

where $(\Delta u_i)_n$ is the elastic opening at time $t_n$; $\delta \Delta u_i^{ve} = (\Delta u_i^{ve})_{n+1} - (\Delta u_i^{ve})_n$ is the incremental interface opening at current time step.

The total cohesion force at time $t_{n+1}$ is finally written as:

$$ c_i^{ve} = c_i + c_i^v = \frac{c_i^{max}(1-D_i)}{D_i} \frac{\Delta u_i^{ve}}{\delta_i^c} + \beta \frac{c_i^{max}(1-D_i)}{D_i \delta_i^c} \left( e^{-\frac{\Delta t}{\lambda_i}} (\Delta u_i^{ve})_n + \frac{1-e^{-\frac{\Delta t}{\lambda_i}}}{\frac{\Delta t}{\lambda_i}} \Delta u_i^{ve} \right) \qquad (22) $$

As with Eq. (14), for the normal part of interface model with negative relative displacement, a penalisation term should be added into Eq. (22).
331 3.3.3. Solid grains

332 The behaviours of quartz, calcite, and pyrite grains are assumed as linear elastic, and the clay aggregates are
333 assumed as viscoplastic. Under the assumption of small strains, the strain tensor of clay aggregates at
334 microscale can be decomposed into its elastic and viscoplastic parts, denoted by superscripts e and $^{vp}$,
335 respectively:

$$
\dot{\epsilon}_{ij} = \dot{\epsilon}_{ij}^{e} + \dot{\epsilon}_{ij}^{vp} \qquad (23)
$$

336 The viscoplastic model for clay aggregates is improved based on Perzyna's overstress concept
337 (Perzyna, 1966) and Lemaitre's creep model (Lemaitre and Chaboche, 1990) by considering creep
338 volumetric strain and a non-associated viscoplastic flow rule. Precisely, the viscoplastic flow rule writes:

$$
\dot{\epsilon}_{ij}^{vp} = \dot{\lambda}^{vp} \frac{\partial G_{vp}}{\partial \sigma_{ij}} \qquad (24)
$$

339 where $\lambda^{vp} \ge 0$ is the viscoplastic multiplier and $G_{vp}$ is the non-associated viscoplastic potential defining
340 the direction of viscoplastic strain rate.

341 Inspired by Lemaitre's model, the viscoplastic multiplier $\dot{\lambda}^{vp}$ is given by (Sun et al., 2023b):

$$
\dot{\lambda}^{vp} = \frac{1}{\eta} \left\langle \frac{F_{vp}}{\sigma_r} \right\rangle^n e^{-k\gamma^{vp}} \quad (25)
$$

342 where $F_{vp}$ is the viscoplastic yield function, $\langle \rangle$ stands for Macaulay's brackets with $\langle x \rangle = \max(0, x)$,
343 while $\eta$, $n$, $k$, and $\sigma_r$ are model parameters. The parameters $n$ (should satisfy $n \ge 1$) and $k \ge 0$
344 represent the effect of stress intensity and strain hardening on the creep strain rate, respectively. The
345 parameter $\eta$ represents the rock viscosity and $\sigma_r$ is a reference stress (1 MPa). The equivalent deviatoric
346 viscoplastic strain $\gamma^{vp}$ is chosen as the hardening parameter, written as $\dot{\gamma}^{vp} = \sqrt{2/3 \dot{e}_{ij}^{vp} \dot{e}_{ij}^{vp}} = \dot{\lambda}^{vp}$, with
347 $e_{ij}^{vp}$ denotes the viscoplastic part of total deviatoric strain tensor $e_{ij}$ which is expressed as $e_{ij} = \varepsilon_{ij} - \varepsilon_{kk}/3 \delta_{ij}$.

349 For the viscoplastic yield function, we suggest the following form:

$$F_{vp} = q + \alpha^{vp} p - \sigma_s \quad (26)$$

350 where $\alpha^{vp} \ge 0$ is a model parameter quantifying the mean stress level sensitivity of the creep behaviour, 351 and $\sigma_s$ is the creep threshold from which viscoplastic strains start to develop. For the viscoplastic potential 352 $G_{vp}$, it takes the following form:

$$G_{vp} = q + \beta^{vp} p \quad (27)$$

353 where $\beta^{vp} > 0$ corresponds to viscous dilatancy of rock and $\beta^{vp} < 0$ corresponds to viscous contraction. 354 Injecting Eqs. (25) – (27) into (24), the viscoplastic strain rate is finally expressed as follows:

$$\dot{\varepsilon}_{ij}^{vp} = \frac{1}{\eta} \left( \frac{q + \alpha^{vp} p - \sigma_s}{\sigma_r} \right)^n e^{-k \gamma^{vp}} \left( \frac{3s_{ij}}{2q} + \frac{\beta^{vp}}{2} \delta_{ij} \right) \quad (28)$$

355 The stresses and tangent operator for elasto-viscoplastic model are updated in each time step at microscale 356 using Newton-Raphson incremental iterative scheme.

357 **3.3.4. Channel flow model**

358 Several assumptions are adopted for microscale hydraulic problems: steady state conditions due to 359 separation of scales, fluid incompressibility, and laminar flow. The water density $\rho_w$ depends on the 360 macroscale water pressure $p_w$:

$$\rho_w = \rho_{w0} \exp\left(\frac{p_w}{K_w}\right) \quad (29)$$

361 where $\rho_{w0}$ (= 1000 kg/m³) is the initial water density at zero fluid pressure; $K_w$ (= 2.2 * 10³ MPa) is 362 the bulk modulus of water.
Considering an interface with two endpoints $s^i$ and $s^j$, the fluid mass flux $\omega$ [kg.m$^{-1}$.s$^{-1}$] of the channel is given by:

$$ \omega = -\phi (p_w(s^j) - p_w(s^i)) \qquad (30) $$

where $p_w(s^i)$ and $p_w(s^j)$ are the fluid pressure at the extremity of each channel, $\phi$ [s] is a hydraulic transmissivity term, written as:

$$ \phi = \frac{\rho_w}{\mu_w} \left( \int_{s^j}^{s^i} \frac{1}{\kappa} ds \right)^{-1} ; \quad \kappa = \frac{\Delta u_h^3}{12} \qquad (31) $$

The variable $\kappa(s)$ [$m^3$] in above equation represents the geometric transmissivity function of the channel with the $s$ coordinate along this channel, formulated with classic cubic relation between the hydraulic equivalent opening of interface $\Delta u_h$, and dynamic viscosity of water $\mu_w$ ($\approx 10^{-9}$ MPa·s). The hydraulic opening $\Delta u_h$ is defined as a piecewise function of mechanical normal opening $\Delta u_n$ (van den Eijnden et al., 2017):

$$ \Delta u_h = \max(\Delta u_h^{\min}, \Delta u_h^{\min} - \Delta u_n^{\min} + \Delta u_n) \qquad (32) $$

where $\Delta u_h^{\min}$ and $\Delta u_n^{\min}$ are introduced to control the initial and minimum permeability of the REA and avoid negative hydraulic conductivity.

## 4. Double-scale numerical simulation at laboratory scale

In this section, the viscous behaviour of clay rock under constant deviatoric stress has been studied at macroscale, on pluricentimetric clay rock specimen. The hydraulic opening parameters between mineral grains are $\Delta u_h^{\min} = 1 \times 10^{-4}$ mm and $\Delta u_n^{\min} = -1 \times 10^{-4}$ mm in Eq. (32), leading to initial major and minor intrinsic permeabilities respectively $5.5 \times 10^{-18}$ m$^2$ and $3.0 \times 10^{-18}$ m$^2$. These large permeabilities allow to neglect the influence of negative pore water pressure. Firstly, the results of macroscale and mesoscale responses are compared to demonstrate the consistency between the models at
two scales (Section 4.1). This comparison is performed for homogeneous macroscale deformation for which identical mesoscale and macroscale behaviours are expected. Secondly, the relation of the shear deformation between meso- and macro- scales is investigated (Section 4.2). Lastly, a complete three-stage creep of clay rock is well reproduced under high deviatoric stress levels and the creep-induced failure process is discussed (Section 4.3). The influence of rock anisotropy on creep behaviour is also studied.

Concerning the clay rock viscosity, at microscale, the viscosity of the clay matrix is implemented in the clay aggregates or in the clay aggregate contacts, or in both. At mesoscale, REAs with a representative length $L_{REA} \approx 100$ µm are used and creep test are reproduced on macroscale samples. It has to be reminded that the mesoscale behaviour is upscaled at macroscale (on each macroscale material point), by scale transition and computational homogenisation, in a double-scale finite element squared (FE²) framework (Section 3). Therefore, in this section, we go beyond the mesoscale REA (i.e. material point) calculations (thoroughly conducted in Sun et al. (2023b)) and conduct 2D plane strain simulations of macroscale BVP using the double-scale numerical model.

The elastic, plastic, and viscous parameters of solid mineral constituents and their damageable contacts have been calibrated using REAs with a representative length $L_{REA} \approx 100$ µm (van den Eijnden et al., 2017; Pardoen et al., 2020; Sun et al., 2023b). They are listed in **Table 2** and **Table 3**, and are used in the following numerical analyses unless otherwise specified.

<table><caption>**Table 2.** Microscale elastic parameters of solid mineral constituents and their damageable contacts.</caption><thead><tr><th>Materials</th><th><i>E</i> (GPa)</th><th><i>ν</i> (-)</th></tr></thead><tbody><tr><td>Tectosilicates (quartz)</td><td>95</td><td>0.074</td></tr><tr><td>Carbonates (calcite)</td><td>84</td><td>0.317</td></tr><tr><td>Heavy minerals (pyrite)</td><td>305</td><td>0.154</td></tr><tr><td>Clay matrix</td><td>2.3</td><td>0.110</td></tr></tbody><thead><tr><th>Interfaces</th><th>δ<sup>c</sup><sub>t/n</sub> (-)</th><th>D<sup>0</sup><sub>t/n</sub> (-)</th><th>c<sup>max</sup><sub>t</sub>(MPa)</th><th>c<sup>max</sup><sub>n</sub>(MPa)</th></tr></thead></table>
<table><caption>**Table 3.** Microscale parameters of viscoplastic clay aggregates and viscoelastic clay aggregate contacts.</caption><tr><td>0.1</td><td>0.001</td><td>2.5</td><td>1.0</td></tr></table>

<table><tr><td>Viscoplastic parameters of clay aggregates</td><td>η (s)</td><td>α<sup>vp</sup></td><td>β<sup>vp</sup></td><td>σ<sub>s</sub> (MPa)</td><td>n</td><td>k</td></tr><tr><td></td><td>9.8*10<sup>12</sup></td><td>0.1</td><td>0</td><td>3</td><td>4.62</td><td>260</td></tr><tr><td>Viscoelastic parameters of clay aggregate contacts</td><td>δ<sub>t/n</sub><sup>c</sup> (-)</td><td>D<sub>t/n</sub><sup>0</sup> (-)</td><td>c<sub>t</sub><sup>max</sup>(MPa)</td><td>c<sub>n</sub><sup>max</sup>(MPa)</td><td>μ<sub>t/n</sub>(MPa s)</td><td>β</td></tr><tr><td></td><td>0.1</td><td>0.001</td><td>0.025</td><td>0.01</td><td>10<sup>7</sup></td><td>100</td></tr></table>

**4.1. Creep behaviour at macroscopic and mesoscopic scales**

**Fig. 8** (a) and (b) shows the schematic representation of the macroscale numerical model in Lagamine finite element code and the mesostructure with 250 numerical cells ($L_{REA} \approx 100$ µm) assigned to each integration point of the macroscale finite elements. A clays rock sample having a width of 25 mm and a height of 50 mm, and thus a slenderness ratio of 2, is considered under biaxial loading creep conditions. The sample is discretised using 10×20 quadrilateral 9-node elements. The tangential displacements are allowed along the top and bottom surfaces of the sample (i.e. smooth boundaries), the vertical displacement of the bottom surface is blocked (i.e. smooth and rigid boundary), and the displacement of the central node of the bottom surface is blocked in both directions to avoid rigid body horizontal translation. Considering the case of saturated rock, drainage is applied on the top and bottom surfaces of the sample. The biaxial creep test is divided in three phases. An isotropic confining pressure of $\sigma_{11}=12$ MPa is firstly applied on the sample external surfaces. Secondly, the vertical loading $\sigma_{22}$ on the sample
416 upper surface is progressively increased to reach the desired values of deviatoric stress $q = \sigma_{22} - \sigma_{11}$
417 and of the ratio $q/q_{max}$, in which $q_{max} = 34.9$ MPa is the statistical mean value of the maximal
418 deviatoric stress from macroscale triaxial creep experiments under 12 MPa of confining pressure
419 (Armand et al., 2017). Then, a creep phase is performed under constant $\sigma_{11}$ and $q$, for $q/q_{max} = 50\%$,
420 75\% and 90\%. At the mesoscale, similar creep tests have been conducted on the REA and used for
421 comparison.

422 Comparisons between the macroscale (double-scale) and the mesoscale creep response are shown in
423 **Fig. 8** (c) and **Fig. 8** (d), in which $\epsilon_{22}$ is the total strain of top surface of rock sample obtained from
424 multiscale numerical simulation, and $\epsilon_{22}^M$ is the total strain calculated from one-point calculation. The
425 viscosity is considered both in clay aggregates and clay aggregate contacts. From **Fig. 8** (c), a good
426 agreement can be observed between the macro- and meso-scale curves for each deviatoric stress ratio
427 $q/q_{max}$. From **Fig. 8** (d), the proposed model is able to well reproduce the progressive decrease of the
428 vertical stress due to relaxation process, and a good agreement can also be found between two curves.

<figure><img src="image_9.png" /></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Double-scale biaxial compression test: (a) schematic representation of the macroscopic configuration and boundary conditions; (b) the mesostructure assigned to each integration point of macroscale finite elements; (c) comparisons of creep simulation at two scales; (d) comparisons of relaxation simulation at two scales.</figcaption></figure>

**4.2. Shear orientation tendency at macroscopic and mesoscopic scales**

Following the analyses in Pardoen et al. (2020), it has been shown that, under deviatoric compression condition, the relative movements between mineral grains (at clay-clay aggregate, clay-inclusion, and inclusion-inclusion contacts) occur mainly by (tangential) slipping along their interfaces, resulting in shear deformations inside the mesostructure. Therefore, an overall mesoscale shear strain of the REA can develop under loading, depending on the loading type (e.g. deviatoric) and on the mesostructure representativeness of the clay rock mechanical behaviour (related to the REA size and the mineral spatial arrangement). Under deviatoric loading and creep, with stress applied perpendicularly to the REA external faces, it is expected that overall strains develop in the principal loading direction (i.e. vertical and lateral). This is the case when the mesostructure is large enough to be representative of the material behaviour. When the considered material mesostructure is smaller and less representative of the material behaviour, mesoscale shear strain $\epsilon_{12}^M$ can develop. **Fig. 9** (a) shows a typical mode of mesostructural
445 shear deformation of the REA. The case $\varepsilon_{12}^M < 0$ corresponds to the overall shear deformation in
446 counterclockwise direction and $\varepsilon_{12}^M > 0$ corresponds to the clockwise direction.

447 To investigate the relationship of shear deformation and shear orientation tendency between macro-
448 and meso- scales during deviatoric creep, we analyse the material response curves of $\varepsilon_{12}^M$ time evolution
449 of several mesostructures. Three different heterogeneous REAs having different sizes have been
450 considered: a representative one of 100×100 µm (the reference mesostructure REA250_1 composed of
451 250 numerical cells, of **Fig. 8** (b)) and two less representative of 50×50 µm (composed of 50 numerical
cells).

452
453 Creep simulations have been conducted, under constant confining pressure of $\sigma_{11}(\sigma_{11}^M) = 12$ MPa and
454 at a deviatoric stress level of $q(q^M)/q_{max} = 75\%$. This confining pressure corresponds to the amplitude
455 of the *in situ* stress at the median depth of the COx formation. The results of mesoscale tests on REAs and of
456 macroscale tests on laboratory samples in terms of creep strain evolution (creep curves) are shown in **Fig. 9** (b). One can observe that the large representative mesostructure exhibits a very small development of
457 shear deformation (both in short- and long-terms) while the small less-representative mesostructures
458 deform in shearing due to their lack of representativity of the material long-term creep behaviour. At the
459 macroscale, similar creep tests have been conducted on clay rock samples and used for comparison. The
460 macroscale creep strain evolutions are shown in **Fig. 9** (b) and have a good consistency with the
461 mesoscale creep curves. The deformations of the samples after 100 days of creep are shown in **Fig. 9** (c).
462
463 From the results of **Fig. 9** (a-c), one can clearly observe that a small bending of the sample is induced by
464 the macroscale shear deformation during the creep test, when considering small REAs. This phenomenon
465 is more obvious under large shear strain. On the other hand, almost no macroscale sample shearing or
466 bending is observed when considering large representative REAs, which is representative of the material
467 macroscale creep behaviour under the considered loading conditions. Moreover, the shear orientation 468 tendency at the two scales is also consistent (Fig. 9 (b)), which means that the shear deformation at 469 macroscale is correctly related to the mesoscale. These results show the feasibility of the meso-macro 470 double-scale FE² computation considering time-dependent material behaviour.

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

(a)

$L_{REA} \approx 30 \ \mu m$ ($\epsilon_{12}^M > 0$)

$L_{REA} \approx 100 \ \mu m$ ($\epsilon_{12}^M > 0$)

(b)

$L_{REA} \approx 30 \ \mu m$ ($\epsilon_{12}^M < 0$)

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

(c)

471 Fig. 9. Shear deformation at two scales: (a) Schematic representation of the shear deformation of the 472 mesostructure (REA); (b) Comparison of time evolution of shear strain $\epsilon_{12}^M(\epsilon_{12})$ of different REAs at 473 two scales; (c) deformation diagrams of rock samples at 100 days (black solid line: initial configuration; 474 red dash line: deformed configuration).

475 **4.3. Creep failure process**

476 The creep strain rate is defined as follows to better describe the different stages of creep process:

$$ \dot{\epsilon}_{ij}^{vp} = \frac{\epsilon_{ij; t_2}^{vp} - \epsilon_{ij; t_1}^{vp}}{t_2 - t_1} = \frac{\Delta\epsilon_{ij}^{vp}}{\Delta t} \quad (33) $$
where $\dot{\varepsilon}_{ij}^{vp}$ is the creep strain rate (either vertical $\dot{\varepsilon}_{22}^{vp}$ or lateral $\dot{\varepsilon}_{11}^{vp}$); $\varepsilon_{ij}^{vp}$, $t_1$, and $\varepsilon_{ij}^{vp}$, $t_2$ are the creep strains at time $t_1$ and $t_2$, respectively, and $\Delta\varepsilon_{ij}^{vp}$ is the creep strain increment.

Considering the double-scale modelling time consumption, a smaller mesostructure with dimension $L_{REA} \approx 50$ µm (with 50 cells, labelled REA50_1) is used instead of the larger mesostructure of **Fig. 8** (b).
This mesostructure is used to reproduce the complete creep behaviour (all stages) under a high constant deviatoric stress. First, a biaxial compression test is simulated to obtain its peak strength under confining pressure of $\sigma_{11} = 12$ MPa at laboratory scale. From the results of deviatoric stress – vertical strain curves, the maximum deviatoric stress that the material can sustain is $q_{max}^{num} = 38.9$ MPa is obtained. Double-scale numerical simulations of macroscale biaxial creep tests are then carried out considering three cases under a constant confining pressure of $\sigma_{11} = 12$ MPa and a constant deviatoric stress of $q = 80\% \cdot q_{max}^{num} = 31$ MPa:

* case 1: viscosity is considered only in the clay aggregates.

* case 2: viscosity is considered only at the contact between clay aggregates.

* case 3: viscosity is considered both in the clay aggregates and at their contacts.

The numerical simulation results until creep failure (maximal creep of the material, see Section 2 and **Fig. 3**) for each viscous case are shown in **Fig. 10** (a)-(c). The results exhibit the clay rock overall creep behaviour by the evolutions of total vertical strain $\varepsilon_{22}$, total lateral strain $\varepsilon_{11}$, and creep rates (Eq. (33)) with time. First, the instantaneous total strains under isotropic and deviatoric compression loading are of about $\varepsilon_{22} = 1.4\%$ vertically and $\varepsilon_{11} = -1.0\%$ laterally. Slightly lower short-term strains (0.1% lower) are obtained if the viscosity is considered only at the contact between clay aggregates (viscous case 2, **Fig. 10** (b)). Then, the long-term time-dependant strain evolution is visible. From **Fig. 10** (a)-(c), a clear three-stage creep process can be observed in both two directions for all the three cases of viscous material behaviour.
During the primary and secondary creep, the magnitude order of the vertical and lateral creep strain rates of both viscous cases 1 and 3 are between 10⁻⁷ and 10⁻⁸ s⁻¹ and 10⁻⁹ s⁻¹ for the viscous case 2. More precisely, the creep rates are: the lower (~1×10⁻⁹ s⁻¹) if viscosity is considered only at the contact between clay aggregates (duration of ~740h to reach the creep failure), slightly larger (~0.5×10⁻⁸ s⁻¹) if viscosity is considered only in the clay aggregates (duration of ~45h to reach the creep failure), and larger (1×10⁻⁷ s⁻¹) if viscosity is considered both in the clay aggregates and at their contacts (duration of ~24h to reach the creep failure). This indicates that considering both viscosities simultaneously accelerates the material creep process. During the tertiary creep stage, the creep rates increase rapidly until reaching a maximum value of about 10⁻⁴ s⁻¹. This accelerated third stage creep leads to the material creep-induced failure. It can be found that the magnitude order of the creep rates in the first two creep stages are similar between numerical simulations and experiments of Liu et al. (2018) and Zhang et al. (2019).

The final total strains at creep failure for all viscous cases are all around $\varepsilon_{22} = 2.2\%$ in the vertical direction and around $\varepsilon_{11} = -2.0\%$ in the lateral direction. One can observe that the generated creep strains are smaller in the primary and second creep if the viscosity is considered only at the contact between clay aggregates (in the viscous case 2, **Fig. 10** (b)); however, the creep strains are larger in the tertiary creep compared to the other two viscous cases (cases 1 and 3, **Fig. 10** (a),(c)), especially the strain differences during creep failure between cases 1 and 2. Another possible way to show this comparison is illustrated in **Fig. 10** (d). A new set of parameters are selected for clay aggregates in case 1 to have similar creep strains in the first two phases compared to case 2. The similar results can be observed from this comparison. When the total vertical strain in case 1 reaches to the strain (1.18%) where the creep failure begins in case 2, it continues to evolve until it reaches about 2.07% (closer to 2% in **Fig. 10** (a)) where creep failure begins. Therefore, the viscosity of the clay aggregate contacts (i.e. the time-dependant viscous sliding between
rigid clay aggregates) seems to have a slightly lower effect on primary and secondary creep stages than the viscosity of the clay aggregates (i.e. the time-dependant overall creep of the clay matrix), but an important contribution to the creep-induced failure process of the clay rock.

Comparing the viscous cases 1 and 2 (Fig. 10 (a)-(b)) to case 3 (Fig. 10 (c)) indicates that the final total vertical and lateral deformations at failure are very close, which indicates that the creep deformation caused by the viscosity in the clay aggregates or at their contacts is not a simple superposition. As long as the micro-damage at the mineral contacts (interfaces) accumulates to a certain extent, the tertiary creep will occur until reaching failure. The similar observations can also be found at mesoscale (Sun et al., 2023b). However, considering both viscosities simultaneously significantly accelerate the material creep process and reduces the time needed to reach the creep failure.

Furthermore, the lateral expansive total strain is lower than the vertical compressive total strain $|ε_{11}| < |ε_{22}|$ under the considered deviatoric loading and creep deformation process, for all creep stages until creep failure. This indicates a global (total) contractive volumetric behaviour. However, the lateral expansive creep strain rate is larger than the vertical compressive creep strain rate $|ε_{11}^{vp}| > |ε_{22}^{vp}|$, indicating an incremental (rate) dilative volumetric behaviour during the creep process. This is different from the results of laboratory observations (Armand et al., 2017; Liu et al., 2018) in which the axial compressive creep strain rates are greater than the lateral expansive ones in all tests, leading to an important volumetric compaction under creeping. The reason for this difference is that the current mesoscale model overestimates the material dilation, which is related to the displacement of the solid mineral constituents and to the opening of the interfaces between them. Since the rearrangement of the solid constituents is not considered in the mesoscale model, this opening at solid constituent contacts occurs regardless of the deformation state. The creep results shown in Fig. 10 (a) and Fig. 10 (c) are very
543

544

545

546

similar except for the time to reach creep failure. Moreover, we only have calibrated viscosity parameters for clay aggregates from mesoscale study (Sun et al., 2023b). Therefore, the viscosity is only considered in the clay aggregates for the rest of the numerical simulations if not specifically stated.

<figure><img src="image_16.png" /><figcaption>(a) Viscosity only of the clay aggregates.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(b) Viscosity only of the contacts between clay aggregates.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(c) Viscosity of both clay aggregates and at their contacts.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(d) Comparison of vertical total strain when the viscosity is considered in clay aggregates or in their contacts.</figcaption></figure>

547

548

**Fig. 10.** Temporal evolution of total strains and creep rates under $\sigma_{11} = 12$ MPa and $q = 80\%$. $q_{max}^{num} = 31$ MPa.

549

550

**Fig. 11** (a) shows the comparisons of the evolution of total vertical strain and creep rate at laboratory (macro $\varepsilon_{22}$ and $\dot{\varepsilon}_{22}^{vp}$) and REA (meso $\varepsilon_{22}^{M}$ and $\dot{\varepsilon}_{22}^{vp,M}$) scales, under constant isotropic confining
pressure of $\sigma_{11}=12$ MPa and deviatoric stress of $q = 31$ MPa. Different from the results shown in **Fig. 8** (c) without the appearance of tertiary creep, the two curves of creep strain, at the two different scales (i.e. mesoscale and macroscale), no longer maintain consistency at larger times at laboratory scale. The vertical creep strain rate at laboratory scale $\dot{\varepsilon}_{22}^{vp}$ shows a clear three-stage variation, whereas the vertical creep strain rate at mesoscale $\dot{\varepsilon}_{22}^{vp,M}$ decreases all the time due to the effect of viscoplastic hardening. The difference between the two curves are gradually obvious after 31h of creep, from the point $P_c$ in **Fig. 11** (a), at which the macroscale creep rate $\dot{\varepsilon}_{22}^{vp}$ begins to increase. Later during the tertiary material creep, the fields of Von Mises' equivalent deviatoric strain $\varepsilon_{eq} = \sqrt{2e_{ij}e_{ij}/3}$ close to the start of the creep failure (point $P_B$) and at creep failure (point $P_A$, end of creep) are shown in **Fig. 11** (b) and (c). At macroscale, deviatoric deformation concentrates in the macroscale sample on the left side of the top surface in this simulation and it is one order of magnitude larger than that in the other locations of the sample. To enlighten the relation between the developments of creep deformation and creep-induced failure across scales, the damage patterns in three deformed mesostructures at different locations in the clay rock sample are shown in **Fig. 11** (c). In **Fig. 11** (c), $D_t^{max}$ and $D_n^{max}$ represent the maximum damage parameter (**Fig. 6**) at tangential and normal directions respectively, and $P(D_t > D_n)$ is the proportion of interfaces in the partially damage state which damage in the tangential direction is greater than that in the opening direction. In the zones where macroscale deviatoric strain develops, it has been observed, at meso- and micro-scales, that damage mainly in shearing develops at the interfaces between mineral grains. It develops mainly and rapidly in the tertiary stage of creep and eventually leads to creep-induced failure (**Fig. 11** (c)). A creep-induced microcrack, represented by full decohesion between mineral grains or in the clay matrix (black parts), develops through the entire mesostructure. However, outside the zone of intense deviatoric
572

strain at macroscale, the material mesostructure is almost not deformed neither damaged. This enlighten a

573

clear relation between creep-induced failure across scales.

<figure><img src="image_20.png" /><figcaption>**Fig. 11.** Creep behaviours of macroscale laboratory sample and of the mesostructure REA50_1 ($L_{REA} \approx 50$ µm), under $\sigma_{11}(\sigma_{11}^M) = 12$ MPa and $q(q^M) = 80\% \cdot q_{max}^{num} = 31$ MPa: (a) comparisons of creep at two scales; macroscale Von Mises' equivalent deviatoric strain field in clay rock sample (b) during tertiary creep (at point $P_B$) and (c) at creep-induced failure (at point $P_A$) with deformed mesostructures.</figcaption></figure>

574

575

576

577

578

The influence of the mesoscopic heterogeneity of the clay rock, from its mineral spatial variability at small scale, on the creep-induced failure at macroscale is studied hereafter. This is done by using another mesostructure (of dimension $L_{REA} \approx 50$ µm and labelled REA50_2) as shown in **Fig. 12**. The creep

579

580
deformation development is studied under a constant confining pressure of $\sigma_{11} = 12$ MPa and deviatoric stress of $q = 30.5$ MPa, corresponding to 85% (stress ratio of 85% is used to make tertiary creep appear) of the mesostructure peak strength $q = 85\% \cdot q_{max}^{num}$. Comparing the creep strain curves of **Fig. 11** (a) and **Fig. 12** (a), obtained for two different material mesostructures, the creep strain during the tertiary creep for the second mesostructure (in **Fig. 12** (a)) is larger. Correspondingly, differences between the development of creep-induced localised deviatoric strain pattern at macroscale and creep-induced localised crack pattern at mesoscale are observed at creep failure and shown in **Fig. 11** (c) and **Fig. 12** (b). For the second mesostructure, a large macroscale shear band, represented by a localised deviatoric strain zone, develops in the bottom of the sample. In this zone of intense macroscale deviatoric strain, a crack develops through the entire mesostructure. Again, this enlighten a clear relation between creep-induced failure across scales, which is observed for different mesostructures. Shi et al. (2021a, 2021b) have observed similar developments of crack patterns through uniaxial compression creep tests on COx clay rock, using X-ray micro-tomography and image analyses.

<figure><img src="image_21.png" /></figure>
<figure><img src="image_22.png" /><figcaption>**Fig. 12.** Creep behaviours of macroscale laboratory sample and of the mesostructure REA50_2 ($L_{REA} \approx 50$ µm) under $\sigma_{11}(\sigma_{11}^M) = 12$ MPa and $q(q^M) = 85\% \cdot q_{max}^{num} = 30.5$ MPa: (a) comparisons of total vertical strains and creep rates at laboratory scale; (b) macroscale Von Mises' equivalent deviatoric strain field in clay rock sample at creep-induced failure and deformed mesostructures.</figcaption></figure>

From the numerical results shown in **Fig. 11** (a), the microcrack development is not visible in the clay rock just before creep failure (point $P_B$ in **Fig. 11** (a)). However, the micro-damage accumulates rapidly with time which lead to the development of microcracking through the entire mesostructure (**Fig. 11** (b)-(c) and **Fig. 12** (b)). This indicates that the creep strain and its rate is an important factor of long-term engineering safety; therefore, the monitoring of creep strain and rate is essential in engineering practice. Tertiary creep control (mitigation) is still difficult because of the high degree of rock damage involved. However, it is crucial to accurately estimate the remaining life time of deep tunnels drilled in rocks after the onset of the tertiary creep. In the absence of precise data, the tunnel should be automatically deactivated to eliminate safety hazards (Frenelus et al., 2022).
Like many sedimentary rocks, the COx clay rock exhibits significant anisotropic properties due to the preferential orientation of quasi-horizontal sedimentary layers and bedding planes (Zhang et al. 2019). The latter can constitute material weakness planes. The petrogenesis processes, involved in the origin and formation of rocks, induce particular microscale morphology of the mineral grains which leads to larger scale anisotropy of the material. Experimental mineralogical and microstructural measurements have been performed in vertical sections perpendicular to the COx clay rock bedding planes (Robinet et al., 2012). The results have shown that quartz (tectosilicates) and calcite (carbonates) mineral inclusions are elongated and have a preferential horizontal average orientation parallel to the bedding plane direction. Furthermore, the carbonate inclusions are smaller, more elongated, and oriented more horizontally (with a less variable orientation) than quartz inclusions (as shown in **Table 1**). Experimental characterisations have emphasised that the macroscale cross-anisotropic behaviour of the COx clay rock structure is related to the orientation and morphology of its mineral inclusions (Robinet et al., 2012). Therefore, these microscale morphological characteristics (i.e. mineral grain elongation, orientation, and size) are accounted for in the generation of the numerical mesostructures and are considered in the mesoscale numerical model (Pardoen et al., 2020).

To investigate the influence of the clay rock anisotropy on its creep behaviour, the anisotropy of the shear strength (i.e. the maximal deviatoric stress) of the clay rock is firstly studied using the mesostructure labelled REA50_2 submitted to a biaxial deviatoric loading. A rotation angle $\theta$ is defined as the orientation of the mesostructure bedding planes with respect to the horizontal direction of the confining stress $\sigma_{11}$ (Fig. 13 (a)); this angle also corresponds to the orientation of the main vertical $\sigma_{22}$ and deviatoric $q = \sigma_{22}-\sigma_{11}$ loading relative to the normal direction to the bedding planes. For instance, an angle of $\theta = 0°$ corresponds to a deviatoric loading perpendicular (⊥) to the bedding planes and an angle of $\theta = 90°$ corresponds to a deviatoric loading parallel (//) to the bedding planes. Such definition of loading direction
relative to the material structure allows to study the anisotropic mechanical response of the material. The peak strengths of the rock mesostructure under different rotation angles $q_{max}^{num}(\theta)$ are shown in **Fig. 13** (b).

It can be clearly seen that the rock shear strength is related to the loading direction. The shear strengths for inclined deviatoric loadings between $\theta = 30^\circ$ and $\theta = 60^\circ$ are significantly lower than those of the two principal material directions, i.e. the shear strength perpendicular to the bedding planes $q_{max}^{num,\perp}$ for $\theta = 0^\circ$ and the shear strength parallel to the bedding planes $q_{max}^{num,//}$ for $\theta = 90^\circ$. This corresponds to a material shear weakness orientation. It can be explained by the fact that the shear failure process of the rock can be enhanced by localised relative sliding between mineral grains along the preferential structural orientation related to bedding structure, mineral grain morphology and orientation. These numerical results on anisotropic shear strength are consistent with experimental observations reported in clay-rich stratified rocks (e.g. Liu et al. 2015; Togashi et al. 2017; Zhang et al. 2019).

After the analysis of the anisotropic material shear strength, the anisotropy of the clay rock creep behaviour is studied. **Fig. 13** (c) shows the evolution of total vertical strain $\varepsilon_{22}$ and vertical creep strain rate $\dot{\varepsilon}_{22}^{vp}$ with different mesostructure rotation angles (i.e. different deviatoric load orientations relative to the bedding planes) of $\theta = 0^\circ, 30^\circ, 60^\circ$, and $90^\circ$ under constant confining $\sigma_{11} = 12$ MPa and deviatoric $q = 21.1$ MPa stresses (80% of the peak strength with $\theta = 60^\circ$, $q = 80\% \cdot q_{max}^{num}(\theta = 60^\circ)$). The rock mesostructure exhibits a creep behaviour for all loading orientations and the mesostructure subjected to a deviatoric constant loading inclined of $\theta = 60^\circ$ deforms unstably up to rupture. On the contrary, when the mesostructure is subjected to the other constant deviatoric loading orientations with, $\theta = 0^\circ, 30^\circ$, and $90^\circ$, it deforms gradually with time. As for the short-term (instantaneous) failure, a material weakness is observed for the long-term (time-dependant) creep-induced failure for an inclined loading direction at $\theta = 60^\circ$ relatively to the material preferential structure.
Furthermore, the time evolutions of the vertical creep strains for oriented loadings of $\theta = 0°, 30°, 60°$ and $90°$ are relatively parallel which corresponds to similar creep rates, as shown in **Fig. 13** (c). The difference between these three strain-time curves comes mainly from the rapid short-term strain increase during the initial deviatoric loading. The creep strain rate for the inclined loading of $\theta = 60°$ during the first two creep stages (primary and secondary creep) is close to that of the other three loading orientation cases. However, it increases with time during tertiary creep stage until creep-induced rupture. This indicates that the anisotropy of the clay rock shear strength and its oriented weakness have an effect on the short-term deformation but only a limited effect on the long-term creep strain development under the constant deviatoric load. These numerical results are consistent with the experimental observations from Zhang et al. (2019). However, the anisotropy of the clay rock seems to have an influence on creep-induced cracking patterns (Shi et al., 2021a), which is out of scope of present study.

<figure><img src="image_23.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(b)</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Fig. 13.** Influence of structural anisotropy on clay rock creep, under $\sigma_{11} = 12$ MPa and $q = 21.1$ MPa: (a) rotation angle $\theta$ defining the orientation of the deviatoric loading relative to the mesostructure (b) peak deviatoric stress for different orientations of the loading; (c) evolution of vertical total strain and creep rate of the mesostructure under loadings oriented of $\theta = 0^\circ$ (⊥), $30^\circ$, $60^\circ$, and $90^\circ$ (//).</figcaption></figure>

**5. Double-scale numerical simulation at gallery scale**

Assessing gallery construction and operational phases are crucial to ensure the stability and sustainability of underground structures, as those required for the deep repository of nuclear wastes. Consequently, the short-term excavation and the induced damaged zone (EDZ) development as well as the long-term convergence evolutions are studied hereafter with the double-scale FE² approach.

At the main level of MHM URL (about -490 m), the *in situ* compressive principal stresses are estimated to be $\sigma_v = 12.7$ MPa in the vertical direction, $\sigma_H = 1.3$ $\sigma_h = 16.1$ MPa in the horizontal direction corresponding to the major principal stress direction, and $\sigma_h = 12.4$ MPa in the horizontal direction corresponding to the minor principal stress direction. The simulated gallery is denoted as GED (Galerie Expérimentale Deux) in the Andra's URL. It has a radius of 2.6 m and is excavated along the minor principal horizontal stress $\sigma_h$, which corresponds to an anisotropic state of stress in the gallery cross section, $\sigma_x = \sigma_H$ and $\sigma_y = \sigma_v$.
## 5.1. Numerical model

A quarter of the gallery is modelled by assuming symmetry along the x and y axes as shown in Fig. 14 (a). The size of the domain and the distance to the outer boundaries is 50 meters. All computations are performed in a two-dimensional plane strain state. Considering the larger deformation, induced by the excavation and evolving with time, in the rock in the vicinity of the gallery, the mesh has been defined with a finer discretisation close to the gallery wall. The mesh has a total of 8181 nodes and 2080 finite elements. The gravity effect is neglected to consider a homogeneous anisotropic initial stress state. The initial pore water pressure field is uniform and equal to $p_{w0} = 4.7$ MPa everywhere corresponding to the in situ pore water pressure at the main level of the MHM URL. The far-field boundaries are under drained conditions with constant pore water pressure $p_{w,0}$ and constant total stresses $\sigma_{x,0} = \sigma_H$ and $\sigma_{y,0} = \sigma_v$ in the normal direction to the boundaries. The symmetry of the model is established by assuming impermeable boundaries and zero normal displacements along the symmetry axes (x and y axes). It should be noted that for the modelling that involves strain localisation to reproduce the EDZ by shear bands and the second gradient model, the normal derivative of radial displacements should be equal to zero along the axis of symmetry. This is because of the existence of gradient terms in the governing equation Eq. (1).

Concerning the excavation operation, the boundary conditions at the gallery wall is modelled by the mechanical and hydraulic deconfinement curves, using the following relations (see Fig. 14 (b)):

$$ \vec{\sigma}_{gal} = (\vec{\sigma}_0 - \vec{\sigma}_{res})(1 - \lambda_m) + \vec{\sigma}_{res} \quad (34) $$

$$ p_{w,gal} = p_{w,0}(1 - \lambda_w) \quad (35) $$

where $\vec{\sigma}_{gal}$ and $p_{w,gal}$ are respectively the stress vector and pore water pressure applied on the gallery wall, $\vec{\sigma}_0$ and $p_{w,0}$ are their corresponding initial values, $\vec{\sigma}_{res}$ is the applied residual stress vector after
704

excavation, $\lambda_m$ and $\lambda_w$ are two deconfinement rate (unloading parameters) which control the stresses and the pore water pressure on the gallery wall (i.e. control the unloading process). As shown in the deconfinement curves of **Fig. 14** (b), the excavation stage is performed in 28 days with $\lambda_w$ evolving from 0 to 1 to reproduce the 3D effect of the advancement of the gallery drilling. After 14 days, the excavation front comes across the studied 2D gallery section and the pore water pressure therefore rapidly drops to zero (Pardoen and Collin, 2017).

705

706

707

708

709

710

After the excavation phase, the long-term operational phase of the unsupported gallery is numerically reproduced. Considering an unsupported gallery (i.e. no gallery wall support), the total stresses and pore water pressure on the gallery wall after the excavation are respectively equal to atmospheric pressure and zero, that is $\sigma_{x,\text{res}} = \sigma_{y,\text{res}} = 0.1 \text{ MPa}$ and $p_{w,\text{gal}} = 0 \text{ Pa}$. These stress and pore water pressure are kept constant at gallery wall in the long term to highlight time-dependant deformation, convergence, and possible induced damage.

711

712

713

714

715

716

<figure><img src="image_26.png" /><figcaption>(a)</figcaption></figure>

717

718
<figure><img src="image_27.png" /><figcaption>**Fig. 14** Schematic representation of the gallery model: (a) numerical model with boundary conditions and finite element mesh; (b) variation of deconfinement rates $\lambda_m$ and $\lambda_w$; (c) random distributions of two mesostructures of dimension $L_{EA} \approx 30$ µm at macroscale.</figcaption></figure>

To consider the spatial heterogeneity and variability of the clay rock composition, different mesostructures containing different mesostructural characteristics are assigned at each macroscale integration point, within the double-scale FE² framework. Mesostructures of dimension $L_{EA} \approx 30$ µm (with 20 numerical cells) are chosen for gallery simulation due to the computational time consumption. Moreover, considering the use of several mesostructures with different mesostructural characteristics (as morphologies and mineral arrangements) lead to different material behaviours which can increase the representativeness of the COx clay rock behaviour at macroscale. A competition study for random macro-variability of different EAs for the gallery excavation can be found in Mourlas et al. (2023). Therefore, two different mesostructures were chosen to increase the representativeness (by variability) of the clay rock behaviour at large-scale. These two mesostructures are respectively labelled EA20_1 and EA20_2 as shown in **Fig. 14** (c) and are randomly assigned at each integration point of macroscale elements.
734 **5.2. Behaviour of selected mesostructures**

735 The parameters used for the microscale constitutive models are adopted from previous calibrations, from
736 **Table 2** and **Table 3**. The mechanical responses of the considered mesostructures (EA20_1 and EA20_2)
737 under biaxial deviatoric compression test perpendicular to the bedding planes and under a confining
738 pressure of $\sigma_{11}$ = 12 MPa are shown in **Fig. 15** (a). Considering hereunder the mesoscale material
739 behaviour, from simulation results, of each mesostructure imply that the interactions between different
740 EAs at macroscale are neglected. It must be emphasised that the comparison to experimental data on
741 macroscale laboratory samples is meaningful only for homogeneous modes of deformation (before the
742 peak-deviatoric stress) (Pardoen et al., 2020). Moreover, due to the existence of possible
743 non-homogeneous deformations, comparing the mesoscale and macroscale deviatoric stress responses
744 only allows a calibration of the micromechanical properties up to the deviatoric stress peak. In **Fig. 15** (a),
745 one can observe that the variability of the material responses is affected by the material mesostructures
746 (Pardoen et al., 2020). For example, the shear strengths of two EAs are respectively of $q_{max}^{num}$ = 35.7
747 MPa and $q_{max}^{num}$ = 36.8 MPa at the vertical strains $\epsilon_{22}^{M}$ = 2% and $\epsilon_{22}^{M}$ = 2.3%. and the initial stiffnesses
748 estimated by the Young's moduli corresponding to the strain $\epsilon_{22}^{M}$ = 0.2% are respectively of $E_{M}$ = 4325
749 MPa and $E_{M}$ = 4686 MPa.

750 After characterising the short-term mechanical response under deviatoric loading, the long-term creep
751 behaviour of the selected mesostructures is evaluated. The simulation results of biaxial creep test with
752 stress ratios $q^M / q_{max}$ of 90%, 75%, and 50% are shown in **Fig. 15** (b). The results highlight the evolution
753 of the creep strain with time as well as the increase of creep deformation under higher deviatoric stress
754 level $q^M / q_{max}$ in accordance to experimental measurements (Armand et al., 2017) and previous modelling
755

(Sun et al., 2023b). Furthermore, a large variability of the rock creep responses can be observed due to small size of the mesostructures (Sun et al., 2023b).

756

757

<figure><img src="image_28.png" /></figure>

EA20_1

(a)

<figure><img src="image_29.png" /></figure>

EA20_2

<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /></figure>

(b)

**Fig. 15.** Short- and long-term mechanical responses of selected clay rock mesostructures used in gallery simulation under: (a) biaxial deviatoric compression test and (b) biaxial deviatoric creep tests with stress ratios of $q^M / q_{max} = 50\%, 75\%, $ and $90\%$.

758

759

760

761

Concerning the hydraulic behaviour at small scale, the hydraulic opening parameters between mineral grains are $\Delta u_h^{min} = 5 \times 10^{-4}$ mm and $\Delta u_h^{tra} = -5 \times 10^{-4}$ mm in Eq. (32), leading to initial major and minor intrinsic permeabilities respectively $4.3 \times 10^{-19}$ m² and $2.9 \times 10^{-19}$ m². The initial intrinsic permeabilities from the numerical model are set larger than the *in situ* permeability

762

763

764
measurements. This is because the permeability evolution due to mesostructure deformation and damage is underestimated in the current mesoscale FE model (and thus at macroscale in the double-scale FE² framework) due to the modelling of the mesostructure in a 2D model (van den Eijnden et al., 2017). Therefore, the initial permeabilities are increased to mitigate this shortcoming.

## 5.3. Gallery excavation and long-term behaviour

The excavation of the GED gallery is numerically reproduced, followed by the long-term operational phase of the unsupported gallery. The hydro-mechanical material behaviour and parameters have been detailed in the previous section.

### 5.3.1. Diametrical convergence

The numerically computed and experimentally observed (in situ) measurements from Armand et al., 2013; Dizier et al., 2023) horizontal and vertical diametrical convergences of the GED gallery are shown and compared in **Fig. 16**. Firstly, the analysis has been performed without considering the creep deformations. The results are shown in dashed lines in **Fig. 16**. In the short-term during the excavation, the simulation results give a good consistency for the vertical diametral convergence of the gallery walls and a different magnitude between the convergences in vertical and horizontal directions corresponding to convergence anisotropy. However, the vertical convergence and the ratio between vertical and horizontal convergences are underestimated. In the long-term, the convergences remain constant due to the end of consolidation, which is not in agreement with the *in situ* observations. To reproduce the gallery convergence evolution in the long term during the operational phase, the viscosity of the clay matrix (i.e. of the clay aggregates) is considered. One can observe that the material creep has a small influence during the short-term excavation stage since similar convergences are obtained, in comparison to the case without material creep.
Nevertheless, the convergences increase rapidly due to the creep, which shows a good agreement with the *in situ* observations, at least in the vertical direction. This demonstrates the importance of incorporating creep in the modelling of the long-term behaviour of clay-rich rock.

<figure><img src="image_33.png" /><figcaption>**Fig. 16.** Evolutions and comparisons of the numerical and experimental gallery diametrical convergences in the: (a) horizontal and (b) vertical directions. Solid line: considering material creep; dashed line: without considering material creep.</figcaption></figure>

### 5.3.2. Strain localisation

The evolution of the shear strain localisation patterns around the gallery during and after the excavation is shown in **Fig. 17**. It is depicted at four moments: at the passage of the front of excavation at the studied section (14 days), at the end of the excavation (28 days), then in the long term after the excavation at 100 days and 1000 days. The shear damage, shear failure, and EDZ induced by the excavation in the rock surrounding the gallery is reproduced by macroscale shear bands (Pardoen et al., 2015a, 2015b; Pardoen and Collin, 2017; Mourlas et al., 2023). The latter correspond to shear strain localisation and accumulation in restricted areas around the gallery and are considered as a precursor to material shear failure. The results are shown in terms of the Von Mises' equivalent deviatoric strain which corresponds to the global amplitude of the shear strain. One can observe that the strain localisation appears slightly during the excavation and keep developing over time. At the end of excavation (at 28 days), shear bands have already
appeared close to gallery wall. Considering that the vertical stress is smaller than the horizontal stress $\sigma_v$ < $\sigma_H$, the fractured zone develops mostly in the vertical direction, as the EDZ observed *in situ* (Armand et al., 2014). In the long term after the excavation, the continuous development of the shear band is due to the development of the creep strain in the clay matrix at microscale, which engenders creep-induced failure at mesoscale (tertiary creep), and creep-induced shear bands which form around the gallery at macroscale. This relation between creep-induced shear strain and failure across scale is even more detail hereafter and is an important result given by multi-scale modelling of time-dependent behaviour and failure of rocks

810

<figure><img src="image_34.png" /></figure>

<figure><img src="image_35.png" /></figure>

<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /><figcaption>**Fig. 17.** Evolution of macroscale shear band (shear strain localisation) pattern representing the EDZ around the GED gallery parallel to the minor horizontal principal stress $\sigma_h$.</figcaption></figure>

813
814

815

816

817

818

819

820

821

822

823

824

825

826

827

828

829

830

831

832

The current activity of the shear deformation process can be show by the Von Mises' equivalent deviatoric strain rates. These fields are again shown at four moments in **Fig. 18**: when the excavation crosses the studied section (14 days), at the end of the excavation (28 days), and in the long term after the excavation at 100 days and 1000 days. It can be observed that the maximum Von Mises' equivalent deviatoric strain rate increases during the excavation whereas it gradually decreases in the long term after the excavation. This is mainly due to the increasingly significant viscoplastic hardening effect. From Eqs. (24) and (25), the hardening term $\exp(-k\gamma^{vp})$ decreases the viscoplastic strain rate $\dot{\epsilon}_{ij}^{vp}$ with the accumulation of deviatoric viscoplastic strain in time. When a shear band is formed in the vicinity of the gallery, the activity of the deformation process varies significantly at different material locations. The development of the rock shear deformation is mainly concentrated within the shear bands at 100 days; however, after a longer time, the deformation is active mainly in one most significant shear band (see the shear band shown in **Fig. 18**). This indicates that the creep of the rock at large times seems to be mainly influenced by the deformation within one (or maybe few) most significant shear band. These dominant shear bands dominate the long-term behaviour of the EDZ around the gallery.
<figure><img src="image_38.png" /><figcaption>**Fig. 18.** Evolution of macroscale shear band activity pattern in the EDZ around the GED gallery parallel to the minor horizontal principal stress $\sigma_h$.</figcaption></figure>

The evolution of deformed and damaged mesostructures around the gallery at different locations inside and outside macroscale shear bands are illustrated in **Fig. 19**. Three moments are chosen to show the material creep effects on the relation between macro- and meso-scale deformations: the end of excavation (28 days) and in the long term after the excavation at 100 days and 1000 days. The development and evolution of localised micro-damage and meso-crack path with time due to clay aggregate viscosity is clearly observed in **Fig. 19**. Creep-induced shear micro-damage at mineral grain contacts leads to shear meso-failure, with meso-shear-cracking developing through the entire EA at mesoscale, are observed. The latter lead to creep-induced shear macro-strain localisation bands (precursor to macro-shear-fracture) at macroscale. Connected meso-crack paths, developing through the entire EA by micro-damage coalescence in microscale localised zones, are visible in the macro-shear band in the
846

847

848

849

850

851

852

853

854

855

856

857

858

859

860

861

862

863

864

865

long term, after 100 and 1000 days, but not at the end of the gallery excavation at 28 days. The
convergence of the host rock (Fig. 16) and the extent of the shear band zone representing the EDZ (Fig.
19) also increases with time due to the clay matrix viscosity. The EDZ refers to the volume of rock in
which crack formation becomes important, thereby deteriorating its mechanical and hydraulic properties.
Consequently, a particular attention should be paid to the EDZ because the evaluation of the damage zone
extension is critical for the design of underground radioactive waste disposal sites (Pellet et al., 2009).
The latter prevent the leakage and migration of radioactive particles (i.e. radionuclides), in the
underground and towards the biosphere, thanks to an engineer and geological multi-barrier system.

Furthermore, a same EA located at different locations at macroscale is subjected to different stress
conditions and can therefore exhibit different meso-crack paths. This can also be due to the interactions
between several EAs at the macroscale and to their variable behaviour, which cannot be reflected at the
scale of only one EA. Furthermore, the interfaces between soft clay aggregates and stiff mineral
inclusions (see the mineral properties in Table 2) seem more likely to be damaged under the excavation
(deviatoric) unloading condition. For instance, a transverse fracture through the entire EA20_1 first
appears above the gallery wall under intense deviatoric stress condition. Moreover, the microscale
damage state of mineral contact interfaces inside a EA develops with time due to the clay aggregates
viscosity, which can play the role of a mesoscale shear crack attractors and then of a macroscale shear
strain localisation attractors, which consequently engenders a larger development zone of shear bands and
EDZ.
<figure><img src="image_39.png" /></figure>

(a)

<figure><img src="image_40.png" /></figure>

(b)

<figure><img src="image_41.png" /></figure>

(c)

866

867

868

869

870

871
Fig. 19. Evolution of macroscale shear band (shear strain localisation) pattern representing the EDZ around the GED gallery and of deformed mesostructures ($L_{EA} \approx 30$ µm) at different locations at (a) the end of excavation (28 days) and after the excavation at (b) 100 days and (c) 1000 days.

## 6. Conclusions

The creep behaviours of the Callovo-Oxfordian (COx) clay rock at several scales are modelled from small-scale viscous mechanisms of the rock medium using double-scale numerical modelling, in a finite element (FE²) framework. To reproduce the viscous behaviour of clay rock, two microscale mechanisms have been introduced in the clay matrix: the viscoplasticity of clay aggregates and the viscoelasticity of their contacts (Sun et al., 2023b). The viscous behaviour of the clay rock has been studied at macroscale, firstly by the modelling of laboratory creep tests on pluricentimetric clay rock samples under deviatoric loading. A consistency of the strain evolution and of the shear strain orientation exist between the two scales, which implies that the macroscale creep behaviour can be reproduced from its mesostructural behaviour. This shows the feasibility of the double-scale computing including time-dependent viscous behaviour. This allows further possibilities of double-scale numerical modelling of large-scale engineering problems, as the prediction of the behaviour of underground structures constructed in rocks. Relation between creep-induced shear strain and failure across scale is an important result given by multi-scale modelling of time-dependent behaviour and failure of rocks in the present study.

The time-dependent creep behaviour of the clay rock is studied under constant deviatoric loading condition as it corresponds to the stress solicitation in the rock surrounding underground galleries in the long term. From the simulation results of laboratory biaxial deviatoric creep tests, a clear three-stage creep process is reproduced. It includes primary creep stage with strain rate deceleration, secondary creep stage with constant creep strain rate, and tertiary creep stage with accelerating strain rate. Such creep behaviour of the COx clay rock has been observed by laboratory triaxial creep tests of Liu et al. (2018). The viscosity
in the clay aggregate contacts has a slight effect on the primary and secondary creeps, but it has an important contribution to the creep-induced failure of the clay rock. This indicates that the microscale damage at the interface between mineral grains develops mainly and rapidly in the tertiary stage, leading eventually to creep failure. The creep strain and its rate is an important factor of long-term engineering safety; therefore, the monitoring of creep strain and rate is essential in engineering practice. The final total strains that develops in the clay rock sample are similar whether the viscosity is considered in the clay aggregates, or in the clay aggregate contacts, or in both. It indicates that as long as the microscale damage at the interface between mineral grains reaches a certain state, the tertiary creep will occur until inducing failure. During tertiary creep, the creep-induced shear strain localisation zone can appear at macroscale, originated from the initiation, growth, accumulation, propagation and coalescence of microcracks and mesocracks within the material at macroscale. Furthermore, the small-scale heterogeneity of the clay rock and the mineral arrangement lead to different rupture patterns at large scale. These results enlighten a clear relation between creep-induced failure across scales. Finally, the anisotropy of the clay rock has limited effect on the creep strain rate.

Then, the excavation and operational phase of the *in situ* Galerie Expérimentale Deux (GED) have been modelled. The simulation results show that the developed multi-scale viscous model is able to simulate the evolution of the long-term gallery convergences due to the viscoplasticity of the clay aggregates. The modelled vertical convergence of the gallery walls has a good consistency compared to *in situ* measurements; however, the horizontal convergence is overestimated. Moreover, the clay aggregate viscosity leads to long-term development of the shear strain localised zone representing the Excavation Damaged Zone (EDZ) around the gallery with time. Different micro-damage states and meso-crack paths are also observed in the rock in the vicinity of the gallery, depending on the macroscale locations. The
917 modelling has also shown that the creep of the rock at large times seems to be mainly influenced by the deformation within one (or maybe few) most significant shear band. These dominant shear bands dominate the long-term behaviour of the EDZ around the gallery. This has been observed through the analyses of macro-shear bands around the gallery and through meso-shear cracks and micro-shear damage induced by the material creep within the shear bands at large scale. From the numerical results at laboratory and gallery scales, the creep-induced shear failure with meso-shear-cracking can develop through the mesostructure, leading to creep-induced shear macro-strain localisation bands (precursor to macro-shear-fracture) at macroscale.

925 **Credit authorship contribution statement**

926 **Yufeng Sun**: Finite element code development, Methodology, Writing – original draft, Writing – review & editing. **Benoit Pardoen**: Finite element code development, Conceptualisation, Writing – review & editing, Supervision. **Henry Kwai-Kwan Wong**: Conceptualisation, Methodology, Writing – review & editing, Supervision.

930 **Declaration of Competing Interest**

931 The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

933 **Data availability**

934 Data will be made available on request.

935 **Acknowledgements**

936 The first author is supported by the China Scholarship Council (CSC) with a PhD Scholarship (File No. 202006560007) for his research work.
938 References

939 Armand, G., Conil, N., Talandier, J., Seyedi, D.M., 2017. Fundamental aspects of the hydromechanical behaviour of Callovo-Oxfordian claystone: From experimental studies to model calibration and validation. Comput. Geotech. 85, 277-286.
941 https://doi.org/10.1016/j.compgeo.2016.06.003.

942 Armand, G., Leveau, F., Nussbaum, C., de La Vaissiere, R., et al., 2014. Geometry and properties of the excavation-induced fractures at the Meuse/Haute-Marne URL drifts. Rock. Mech. Rock. Eng. 47 (1), 21-41.
944 https://doi.org/10.1007/s00603-012-0339-6.

945 Armand, G., Noiret, A., Zghondi, J., Seyedi, D.M., 2013. Short- and long-term behaviors of drifts in the Callovo-Oxfordian claystone at the Meuse/Haute-Marne Underground Research Laboratory. J. Rock Mech. Geotech. Eng. 5 (3), 221-230.
947 https://doi.org/10.1016/j.jrmge.2013.05.005.

948 Bertrand, F., Buzzi, O., Bésuelle, P., Collin, F., 2020. Hydro-mechanical modelling of multiphase flow in naturally fractured coalbed using a multiscale approach. J. Nat. Gas. Sci. Eng. 78, 103303. https://doi.org/10.1016/j.jngse.2020.103303.

950 Chen, C., Xu, T., Heap, M.J., Baud, P., 2018. Influence of unloading and loading stress cycles on the creep behavior of Darley Dale Sandstone. Int. J. Rock Mech. Min. Sci. 112, 55-63. https://doi.org/10.1016/j.ijrmms.2018.09.002.

952 Collin, F., Chambon, R., Charlier, R., 2006. A finite element method for poro mechanical modelling of geotechnical problems using local second gradient models. Int. J. Numer. Meth. Eng. 65 (11), 1749-1772. https://doi.org/10.1002/nme.1515.

954 Cosenza, P., Prêt, D., Giraud, A., Hedan, S., 2015a. Effect of the local clay distribution on the effective elastic properties of shales. Mech Mater. 84, 55-74. https://doi.org/10.1016/j.mechmat.2015.01.016.

956 Cosenza, P., Prêt, D., Zamora, M., 2015b. Effect of the local clay distribution on the effective electrical conductivity of clay rocks. J. Geophys. Res. Solid. Earth. 120, 145-168. https://doi.org/10.1002/2014JB011429.

958 Desbois, G., Höhne, N., Urai, J.L., Bésuelle, P., Viggiani, G., 2017. Deformation in cemented mudrock (Callovo-Oxfordian Clay) by microcracking, granular flow and phyllosilicate plasticity: insights from triaxial deformation, broad ion beam polishing and scanning electron microscopy. Solid. Earth. 8 (2), 291-305. https://doi.org/10.5194/se-8-291-2017.

961 Desrues, J., Argilaga, A., Caillerie, D., et al., 2019. From discrete to continuum modelling of boundary value problems in geomechanics: An integrated FEM-DEM approach. Int. J. Numer. Anal. Meth. Geomech. 43 (5), 919-955.
963 https://doi.org/10.1002/nag.2914.

964 Dizier, A., Scibetta, M., Armand, G., Zghondi, J., et al., 2023. Stability analysis and long term behaviour of deep tunnels in clay formation. Geol. Soc. Spec. Publ. 536. https://doi.org/10.1144/SP536-2022-86.

966 van den Eijnden A.P, Bésuelle, P., Chambon, R., Collin, F., 2016. A FE² modelling approach to hydromechanical coupling in cracking-induced localization problems. Int. J. Solids. Struct. 97-98, 475-488. https://doi.org/10.1016/j.ijsolstr.2016.07.002.

968 van den Eijnden A.P., Bésuelle, P., Collin, F., Chambon, R., Desrues, J., 2017. Modeling the strain localization around an underground gallery with a hydro-mechanical double scale model; effect of anisotropy. Comput. Geotech. 85, 384-400.
970 https://doi.org/10.1016/j.compgeo.2016.08.006.

971 Fabre, G., Pellet, F., 2006. Creep and time dependent damage in argillaceous rocks. Int. J. Rock Mech. Min. Sci. 94, 55-63. 43 (6), 950-960. https://doi.org/10.1016/j.ijrmms.2006.02.004.

973 Frenelus, W., Peng, H., Zhang, J.Y., 2022. Creep behavior of rocks and its application to the long-terms of deep rock tunnels. Appl. Sci. 12 (17), 8451. https://doi.org/10.3390/app12178451.

975 Frey, J., Chambon, R., Dascalu, C., 2013. A two-scale poromechanical model for cohesive rocks. Acta. Geotech. 8 (2), 107-124.
976 https://doi.org/10.1007/s11440-012-0173-8.
977

Fu, T.F., Xu, T., Heap, M.J., Meredith, P.G., Mitchell, T.M., 2020. Mesoscopic time-dependent behavior of rocks based on three-dimensional discrete element grain-based model. Comput. Geotech. 121, 103472. https://doi.org/10.1016/j.compgeo.2020.103472.

978

979

980

Hill, R., 1965. A self-consistent mechanics of composite materials. J. Mech. Phys. Solids. 13 (4), 213-222. https://doi.org/10.1016/0022-5096(65)90010-4.

981

982

Jung, S., Vu, M.N., Pouya, M., Ghabezloo, S., 2022. Effect of anisotropic creep on the convergence of deep drifts in Callovo-Oxfordian claystone. Comput. Geotech. 152, 105010. https://doi.org/10.1016/j.compgeo.2022.105010.

983

984

Lemaitre, J., Chaboche, J., 1990. Mechanics of Solid Materials. Cambridge University Press, Cambridge.

985

Liu, Z.B., Xie, S.Y., Shao, J.F., Conil, N., 2015. Effects of deviatoric stress and structural anisotropy on compressive creep behavior of a clayey rock. Appl. Clay. Sci. 114, 491-496. https://doi.org/10.1016/j.clay.2015.06.03.

987

Liu, Z.B., Shao, J.F., Xie, S.Y., Conil, N., Zha, W.H., 2018. Effects of relative humidity and mineral compositions on creep deformation and failure of a claystone under compression. Int. J. Rock. Mech. Min. 103, 68-76. https://doi.org/10.1016/j.ijrmms.2018.01.015.

988

989

990

Mandel, J., 1972. Plasticité classique et viscoplasticité, vol CISM Lecture Notes. Springer, New York, pp. 187.

991

Mánica, M.A., Gens, A., Vaunat, J., Armand, J., Vu, M., 2021. Numerical simulation of underground excavations in an indurated clay using non-local regularisation. Part 1: formulation and base case. Géotechnique. 72 (12), 1092-1112. https://doi.org/10.1680/jgeot.20.P.246.

992

993

994

Marinelli, F., van den Eijnden, A.P., Sieffert, Y., Chambon, R., Collin, F., 2016. Modeling of granular solids with computational homogenization: comparison with Biot's theory. Finite Elem. Anal. Des. 119, 45-62. https://doi.org/10.1016/j.finel.2016.05.003.

995

996

997

Mourlas, C., Pardoen, B., Bésuelle, P., 2023. Large-scale failure prediction of clay rock from small-scale damage mechanisms of the rock medium using multiscale modelling. Int. J. Numer. Anal. Methods Geomech. 1-35. https://doi.org/10.1002/nag.3513.

998

999

1000

Ozdemir, I., Brekelmans, W., Geers, M., 2008. Computational homogenization for heat conduction in heterogeneous solids. Int. J. Numer. Methods Eng. 73 (2), 185-204. https://doi.org/10.1002/nme.2068.

1001

1002

Pardoen, B., Bésuelle, P., Dal Pont, S., Cosenza, P., Desrues, J., 2020. Accounting for small-scale heterogeneity and variability of clay rock in homogenised numerical micromechanical response and microcracking. Rock. Mech. Rock. Eng. 53 (6), 2727-2746. https://doi.org/10.1007/s00603-020-02066-7.

1003

1004

1005

Pardoen, B., Collin, F., 2017. Modelling the influence of strain localisation and viscosity on the behaviour of underground drifts drilled in claystone. Comput. Geotech. 85, 351-367. https://doi.org/10.1016/j.compgeo.2016.05.017.

1006

1007

Pardoen, B., Seyedi, D.M., Collin, F. 2015a. Shear banding modelling in cross-anisotropic rocks. Int. J. Solids. Struct. 72, 63-87. https://doi.org/10.1016/j.ijsolstr.2015.07.012.

1008

1009

Pardoen, B., Levasseur, S., Collin, F. 2015b. Using local second gradient model and shear strain localisation to model the excavation damaged zone in unsaturated claystone. Rock. Mech. Rock. Eng. 48, 691-714. https://doi.org/10.1007/s00603-014-0580-2.

1010

1011

1012

Pellet, F., Roosevelt, M., Deleruyelle, F., 2009. On the 3D numerical modelling of the time-dependent development of the damaged zone around underground galleries during and after excavation. Tunn. Undergr. Space Technol. 24 (6), 665-674. https://doi.org/10.1016/j.tust.2009.07.002.

1013

1014

1015

Perzyna, P., 1966. Fundamental Problems in Viscoplasticity. Advances in Applied Mechanics, Vol. 9. Academic Press.

1016

Rahal, S., Sellier, A., Casaux-Ginestet, G., 2017. Poromechanical consolidation and basic creep interactions around tunnel excavation. Int. J. Rock Mech. Min. Sci. 94, 55-63. https://doi.org/10.1016/j.ijrmms.2016.12.009.

1017
1018

1019

1020

1021

1022

1023

1024

1025

1026

1027

1028

1029

1030

1031

1032

1033

1034

1035

1036

1037

1038

1039

1040

1041

1042

1043

1044

1045

1046

1047

1048

1049

1050

1051

1052

1053

1054

1055

1056

1057

1058

Robinet, J.C., 2008. Minéralogie, porosité et diffusion des solutés dans l'argilite du Callovo-Oxfordien de Bure (Meuse/Haute-Marne, France) de l'échelle centimétrique à micrométrique. PhD thesis, Université de Poitiers, Poitiers.

Robinet, J.C., Sardini, P., Coelho, D., Parneix, J.C., Prêt, D., Sammartino, S., Boller, E., Altmann, S., 2012. Effects of mineral distribution at mesoscopic scale on solute diffusion in a clay-rich rock: Example of the Callovo-Oxfordian mudstone (Bure, France). Water Resour. Res. 48(5). https://doi.org/10.1029/2011WR011352.

Seyedi, D.M., Armand, G., Noiret, A., 2017. "Transverse Action" – A model benchmark exercise for numerical analysis of the Callovo-Oxfordian claystone hydromechanical response to excavation operations. Comput. Geotech. 85, 287-305. https://doi.org/10.1016/j.compgeo.2016.08.008.

Shi, H.L., Hosdez, J., Rougelot, T., et al., 2021a. Influences of structural anisotropy and heterogeneity on three-dimensional strain fields and cracking patterns of a clay-rich rock. Acta. Geotech. 16, 2175-2187. https://doi.org/10.1007/s11440-021-01152-6.

Shi, H.L., Hosdez, J., Rougelot, T., et al., 2021b. Analysis of local creep strain field and cracking process in claystone by X-Ray micro-tomography and digital volume correlation. Rock. Mech. Rock. Eng. 54:1937-1952. https://doi.org/10.1007/s00603-021-02375-5.

van der Sluis, O., Schreurs, P.J.G., Brekelmans, W.A.M., Meijer, H.E.H., 2000. Overall behaviour of heterogeneous elastoviscoplastic materials: effect of microstructural modelling. Mech. Mater. 32 (8), 449 - 462. https://doi.org/10.1016/S0167-6636(00)00019-3.

Shahbodagh, B., Mac, T.N., Esgandani, G.A., Khalili, N., 2020. A bounding surface viscoplasticity model for time-dependent behavior of soils including primary and tertiary creep. Int. J. Geomech. 20(9), 04020143. https://doi.org/10.1061/(ASCE)GM.1943-5622.0001744.

Souley, M., Armand, G., Kazmierczak, J.B., 2017. Hydro-elasto-viscoplastic modeling of a drift at the Meuse/Haute-Marne underground research laboratory (URL). Comput. Geotech. 85, 306-320. https://doi.org/10.1016/j.compgeo.2016.12.012.

Sun, Y., Wong, H., Pardoen, B., Deleruyelle, F., Dufour, N., Branque, D., Leo, C. 2021. Analytical study of post-closure behaviour of a deep spherical cavity in a dilatant viscoplastic rock mass. Comput. Geotech. 139 (3), 104419. https://doi.org/10.1016/j.compgeo.2021.104419.

Sun, Y., Wong, H., Pardoen, B., Deleruyelle, F., Dufour, N. 2023a. Long-term hydromechanical behaviour of a deep cavity taking into account a simplified sequence of life stages. Comput. Geotech. 156, 105256. https://doi.org/10.1016/j.compgeo.2023.105256.

Sun, Y., Pardoen, B., van den Eijnden, A.P., Wong, H., 2023b submitted. Modelling the time-dependent mechanical behaviour of clay rocks based on meso- and micro-structural viscous properties. Int. J. Numer. Anal. Methods Geomech.

Terada, K., Hori, M., Kyoya, T., and Kikuchi, N., 2000. Simulation of the multi-scale convergence in computational homogenization approaches. Int. J. Solids. Struct. 37 (16), 2285-2311. https://doi.org/10.1016/S0020-7683(98)00341-2.

Togashi, Y., Kikumoto, M., Tani, K., 2017. An experimental method to determine the elastic properties of transversely isotropic rocks by a single triaxial test. Rock. Mech. Rock. Eng. 50 (1), 1-15. https://doi.org/10.1007/s00603-016-1095-9.

Weng, X.L., Li, H., Hu, J.B., Li, L., 2023. Behavior of saturated remolded loess subjected to coupled change of the magnitude and direction of principal stress. Int. J. Geomech. 23 (1), 04022244. https://doi.org/10.1061/(ASCE)GM.1943-5622.0002612.

Xu, T., Zhou, G.L., Heap, M.J., Yang, S.Q., Konietzky, H., Baud, P., 2018. The modeling of time-dependent deformation and fracturing of brittle rocks under varying confining and pore pressures. Rock. Mech. Rock. Eng. 51, 3241-3263. https://doi.org/10.1007/s00603-018-1491-4.
1059

1060

1061

1062

1063

1064

1065

1066

1067

1068

1069

1070

1071

1072

1073

Yang, L., Sheng, X.C., Zhang, B., 2023. Propagation behavior of hydraulic fractures in shale under triaxial compression considering the influence of sandstone layers. Gas. Sci. Eng. 110, 204895. https://doi.org/10.1016/j.jgsce.2023.204895.

Yuan, Y., Xu, T., Heap, M.J., Meredith, P., Yang, T.H., Zhou, G.L., 2021. A three-dimensional mesoscale model for progressive time-dependent deformation and fracturing of brittle rock with application to slope stability. Comput. Geotech. 135, 104160. https://doi.org/10.1016/j.compgeo.2021.104160.

Zalamea, N., Bésuelle, P., Dal Pont, S., Di Donna, A. 2021. Multiscale modelling of the thermo-hydromechanical behaviour of argillaceous rocks. 16th International Conference on Computational Plasticity. Fundamentals and Applications (COMPLAS), Barcelona.

Zener, C., 1948. Elasticity and Anelasticity of Metals. University of Chicago Press, Chicago.

Zhan, Y., Shao, J.F., Duveau, G., Vu, M.N., Armand, G., 2021. Numerical modeling of deformation and damage around underground excavation by phase-field method with hydromechanical coupling. Comput. Geotech. 138, 104369. https://doi.org/10.1016/j.compgeo.2021.104369.

Zhang, C.L., Armand, G., Conil, N., Laurich, B., 2019. Investigation on anisotropy of mechanical properties of Callovo-Oxfordian claystone. Eng. Geol. 251, 128-145. https://doi.org/10.1016/j.enggeo.2019.02.008.