Modelling a gas injection experiment incorporating embedded fractures and heterogeneous material properties

This study focuses on the modelling of a gas injection experiment to assess the effects of incorporating heterogeneous material properties. The numerical model considers a two-phase flow coupled hydro-mechanical problem, and includes embedded fractures that open with deformation, thereby enhancing permeability. The approach used is integrated in the CODE_BRIGHT software, which allows for the consideration of geomaterials with a spatially correlated heterogeneous field of porosity that follows a normal distribution. This spatial correlation can be either isotropic or anisotropic. A key aspect of this approach is that material properties such as intrinsic permeability, diffusivity or cohesion are defined as a function of porosity. Consequently, these properties also exhibit heterogeneity with spatial correlation and, eventually, anisotropy. The results derived from the numerical model align well with in-situ measurements. The study also includes sensitivity analyses to the variation of critical variables. The calibration of the model has been validated through a similar experiment. The findings indicate that the consideration of heterogeneous material properties can have a significant influence on gas injection problems, particularly when a hydraulic fracture is formed.


Introduction
The work reported in this article is part of the DECOVALEX (DEvelopment of COupled models and their VALidation against Experiments) project, an international collaboration focused on advancing the understanding and modeling of coupled thermo-hydro-mechanicalchemical (THMC) processes in nuclear waste repositories and other geo-engineering applications.Task A of DECOVALEX-2023 specifically investigates the behavior of Callovo-Oxfordian claystone (COx) under fluid pressurization and the resulting fracturing.
Callovo-Oxfordian (COx) claystone, a sedimentary clayey rock from the Paris basin deposited 155 million years ago into a 150-m-thick layer slightly tilting with an inclination of 2 • with respect to horizontal over a layer of Dogger limestone and covered by a layer of Oxfordian limestone, has been selected as a potential host rock for a radioactive waste repository by the French agency for the management of radioactive waste (Andra). 1,2nderstanding gas migration within the COx is imperative for evaluating repository performance and its long-term evolution.If the gas generation rate surpasses the rate of dissolved gas diffusion, pressure build-up may exceed hydraulic pressure and capillary resistance in the surrounding rock, potentially causing damage to the rock mass and jeopardizing the safety function of the host rock.
The focus of this article is on modelling of the PGZ experiment 3,4 conducted at the Meuse/Haute-Marne Underground Research Laboratory (MHM URL), carried out to investigate gas transfer in the COx.
It is widely recognized that materials in the geosphere display heterogeneity, meaning that their properties vary spatially.For instance, parameters such as intrinsic permeability may exhibit spatial variation due to the distribution and properties of constituent minerals. 5This heterogeneity could be significant when modelling gas injection problems, as it could influence both pressure build-up and gas migration pathways.
Geostatistics, 6 particularly semi-variogram analysis and stochastic simulation, can be used to describe spatial variability, which is a fundamental characteristic of many natural materials.The software CODE_BRIGHT, 7 developed at the Department of Civil and Environmental Engineering at the Universitat Politècnica de Catalunya (UPC) and utilized in combination with the pre/post-processor GID, 8 developed by the International Centre for Numerical Methods in Engineering (CIMNE), incorporates an approach for modeling heterogeneous material properties.][14][15][16] The general formulation of CODE_BRIGHT is presented in detail in Olivella 17 and Olivella et al., 7 including the weak forms of the governing equations and explicit definitions of the resulting matrices and vectors.Olivella et al. 7 includes explanations of how variables are specifically treated in the nodes, elements and cells.
While heterogeneity is commonly considered for specific parameters for various problems, such as transmissivity for hydrogeological problems, 18,19 strength for geotechnical problems 20,21 or thermal conductivity for heat flow problems, 22,23 its introduction in gas injection problems within a two-phase flow coupled hydro-mechanical (HM) framework is rare.
Several models have been developed to simulate gas flow through low permeability clay formations, such as Callovo-Oxfordian (COx) claystone, 24,25 Opalinus Clay (OPA) 26 and Boom Clay. 27However, these models do not consider heterogeneity of material properties.Therefore, one of the aims of this article is to assess the impact of considering heterogeneous material properties on the results.
An important aspect of this approach is using a variogram to account for spatial correlation, which has been previously verified by Rodriguez-Dono et al. 28 The validation of this approach itself, through comparison with experimental results, has also been presented by Rodriguez-Dono et al. 28 Further information about the features and capabilities of this approach in modeling heterogeneous geomaterials can be found in Rodriguez-Dono et al. 28 Note that the approach used may be especially relevant when dealing with clay-based geomaterials, although it may have limitations when dealing with low porosity geomaterials or in simulating large-scale structures.

Theoretical background
The spherical semi-variogram has been selected as the foundation for the implementation in CODE_BRIGHT, 7 as it is perhaps the most frequently used semi-variogram model. 6,29In the spherical semi-variogram, C 0 represents the nugget, which is attributed to both measurement errors and micro-variabilities of the mineralization.The term a denotes the range, implying that any data value will correlate with any other value within a radius a (as depicted in Fig. 1).
The phenomenon under study is termed "anisotropic" when the modulus of the examined property h varies in different directions. 29In such cases, the directional graph of the range a forms an ellipsoid in   30 and modified from Yven et al. 31 A. Rodriguez-Dono et al.
three dimensions (3D).To derive the anisotropic covariance function, the first step involves rotating the coordinate axes by the angles α, β, and γ (Fig. 2), aligning them parallel to the primary axes of the ellipsoid.Subsequently, the axes of the ellipsoid are defined, coinciding with the ranges in the three principal directions.The major range of the ellipsoid is a x (shown in Fig. 2) and the other ranges are obtained by the ratios λ 1 = a y /a x and λ 2 = a z /a x , where a x , a y , and a z are the ranges for the rotated axes x, y, and z, respectively.More details on this approach can be found in Rodriguez-Dono et al. 28 In the CODE_BRIGHT formulation, many parameters can be functions of porosity, which means heterogeneous characteristics are automatically inherited if porosity is considered heterogeneous.In addition, it was discovered that the pore size distribution of some rocks (e.g.COx, a clay rock), approximately follows a normal distribution (Armand et al. 1 ; Fig. 3b).These reasons led to the selection of porosity as the variable used to generate a spatially correlated field following a normal distribution.
A specific module has been developed in CODE_BRIGHT to generate a spatially correlated heterogeneous field of porosity, 28 following a Gaussian random distribution of mean μ and standard deviation σ (Fig. 3a).The various porosity values obtained are subsequently applied to the elements in the finite element mesh.In the practical use of this approach, only values of the distribution between μ − 3σ and μ +3σ are considered to avoid numerical difficulties related to outliers.
Furthermore, the range a corresponding to the variogram model should be at least twice the size of the mesh elements in order to maintain some spatial correlation.If the range is lower than that, no spatial correlation is observed.
The advective flux q f related to fluid motion is governed by the generalized Darcy's law, as shown in Eq. 1, where k represents the intrinsic permeability tensor; k r represents relative permeability; p f is fluid pressure; ρ f and μ f are density and viscosity of the fluid, respectively; and g represents the vector of gravitational forces.
In modelling complex fluid injection problems, intrinsic permeability can be considered as the sum of two terms (Eq.2).The first term (k matrix ) is a function of the porosity ϕof the matrix of the material (Fig. 4).The second term (k fracture ) corresponds to the permeability associated with the aperture b of embedded fractures, according to a cubic law (Eq.3), [32][33][34] where s is the spacing among fractures (Fig. 5).These small micro-fractures or openings may develop in the material leading to preferential flow paths.Thus, this approach permits differentiating between the intrinsic permeability associated with the matrix and with the fractures. 35 Intrinsic permeability of the matrix can be modelled using an exponential law (Fig. 4), as in Eq. 4, where k 0 represents the intrinsic permeability corresponding to a reference porosity ϕ 0 of the matrix, and b k is a calibration coefficient.Note that the exponential law has been found to represent more accurately the variation in intrinsic permeability observed in certain geomaterials such as bentonite. 15However, it requires calibration due to its empirical nature.The aperture b is calculated as shown in Eq. 5, where ε and ε 0 refer to current and initial strains, respectively (with ε > 0 meaning extension), whereas b 0 and b max are the initial and maximum apertures, respectively.
The retention curve is defined by Eq. 6, where P g , P l are the gas and liquid pressures according to Van Genuchten's model; λ c is a shape function; S e is the effective liquid saturation; S l , S rl and S ls are liquid saturation, residual liquid saturation and maximum liquid saturation, respectively; P o is the reference capillary pressure and σ o is the reference surface tension.
Regarding the relative permeability of the gas phase, a generalized power equation (Eq.7) may be used, 36 where S eg is the effective gas saturation; S g , S rg and S gs are gas saturation, residual gas saturation and maximum gas saturation, respectively; A and λ g are model parameters.Note that A is a gas enhancement factor which is attributed to changes in clay material structure.
Diffusion is a process modelled using Fick's law (Eq.8), where i f is the diffusive flux of a component in a fluid, ϕ is porosity, ρ f is the density of the fluid, S f is the saturation degree of the fluid, τ is the coefficient of Fig. 4. Experimental data corresponding to bentonite and adopted models for the intrinsic permeability law, modified from Sánchez et al. 15 Note that the slope of the exponential law can be adjusted while Kozeny's law has no parameter to allow adjustment.(adapted from Olivella et al. 33 ).
tortuosity, D m is the molecular diffusion coefficient and ω f is the mass fraction of the component in the fluid.It can be applied, for example, to dissolved air.Note I is the identity matrix and D eff is the coefficient of effective diffusion.In the case of air dissolved in the liquid phase, the corresponding molecular diffusion coefficient D air m can be calculated as in Eq. 9, where T is temperature in ºC, D air = 1.1 × 10 − 4 m 2 s − 1 and Q= 24530Jmol -1 .R is the ideal gas constant (8.314JK − 1 mol − 1 ).For more details, refer to Olivella. 17 It can be noted that the diffusive flux i f depends on porosity.Therefore, a heterogeneous distribution of porosity will result in a heterogeneous distribution of diffusivity.Tortuosity is usually considered constant and can be used as a calibration parameter.With constant tortuosity, the equivalent diffusivity depends linearly on porosity ϕ.Accordingly, the coefficient of effective diffusion can undergo variations in space and time.
In the study of material mechanical parameters, certain models incorporate internal variables associated with deformation.][39][40][41] These models decompose the total strain into two components: the elastic strain tensor ε e and the inelastic strain tensor ε i , as shown in Eq. 10.The increments in elastic strain can be a combination of elastic volumetric strain ε e v (Eq.11) and elastic deviatoric strain ε e d (Eq.12), where is the shear modulus, and ν is the Poisson's ratio.
In this article, the effective stress is defined (for compression positive) as σ ′ = σ − max , where σ is a normal component of the total stress tensor and P g , P l are gas and liquid pressures, respectively.The viscoplastic model proposed by Song et al. 41 adopts a Drucker-Prager yield function (Eq.13), where q is the equivalent deviatoric stress, p ′ = (σ ′ x +σ ′ y +σ ′ z )/3 is the mean effective stress, c is the cohesion; M and β are parameters that can be given in compression (Eq.14) and extension (Eq.15) conditions, where ϕ ′ is the friction angle.
F < 0 indicates an elastic stress state, while F ≥ 0 signifies viscoplastic behavior, as given by Eq. 16, where t is time, Γ is a fluidity parameter, Φ(F) = F m is a stress function with a constant power m, σ ′ is the effective stress tensor and G f is the flow rule defined in Eq. 17, where α is a non-associativity parameter that ranges from 0 to 1. ) However, determining the accumulated plastic strain in a sample prior to testing remains an unresolved issue.While we can measure strain by testing a sample in the lab, we cannot measure the historical strain that the sample has undergone before testing.This issue underscores the potential value of introducing porosity as a measure of the material's internal structure since porosity can be measured independently of the material's history.In fact, creep laws in which deformation is a function of void ratio have a consistent theoretical justification. 42n terms of crack formation due to desiccation, cohesion c has been considered as a function of porosity, 43 as shown in Eq. 18, where s is suction, ϕ c0 is a reference porosity and n c , a 1 , a 2 are parameters to be calibrated with experimental data for a particular geomaterial.
By incorporating Eq. 18, the viscoplastic model inherits heterogeneous properties from the heterogeneous porosity distribution.Thus, representing cohesion as a function of porosity may be useful for representing fractures or weakness zones while keeping the friction value close to its original value.It is important to note that cohesion would vary both spatially and temporally following this approach.Moreover, experimental results have indicated that an elastic-viscoplastic approach is more adequate than an elastic-plastic approach for reproducing the yielding behavior of geomaterials. 17The viscoplastic approach also facilitates control over the size of localized areas and thus avoids dependence on the mesh used. 38,41,44The elastic-viscoplastic model has been frequently applied to reproduce time-dependent effects of geomaterials. 38

Modelling a short-term gas flow injection experiment
In this section, a short-term gas injection experiment is modelled and a series of sensitivity analyses are conducted: -Sensitivity to heterogeneous properties.
-Sensitivity to the initial gas saturation degree in the interval.
-Sensitivity to the interval volume.
-Sensitivity to the constitutive model chosen: elastic versus plastic; isotropic versus anisotropic.-Sensitivity to the embedded fractures spacing.

Geometry, initial and boundary conditions, and material parameters
Using CODE_BRIGHT, a two-dimensional (2D) coupled two-phase flow hydro-mechanical model is considered to simulate a gas flow injection experiment.The modelled experiment was conducted in a borehole named PGZ1003, drilled in February 2020, oriented parallel to the major horizontal stress.The borehole is equipped with a multiple packer system to monitor water/gas pressure across five isolated intervals.Between intervals, the borehole annulus is filled with a resin that insures long-term stability and very low compressibility.
During gas injection, the injection interval is isolated by two waterinflated packers connected to a high-pressure vessel to stabilize induced pressure fluctuations.The borehole is sub-horizontal, with a length of approximately 35 m.The length of the packers is 1.5 m and the length of the intervals is 1 m.Each interval comprises a borehole annulus, a cylindrical porous sinter with 30% porosity, and a technological void connected to two hydraulic lines (Fig. 6).These lines facilitate pressure measurement and water or gas injection.The borehole diameter measures 76 mm.
The model considers one quarter of a perpendicular section to the borehole, located at the mid-section of the injection interval #4 (Fig. 7).This section is a square of dimensions 10 m × 10 m, with an additional 1 m considered in the z direction (perpendicular to the modelled section).The initial radius of the modelled borehole is 0.038 m, and the internal radius of the injection interval is 0.033 m.The model considers two different materials: the Callovo-Oxfordian claystone (COx) where the borehole is excavated, and the interval where gas is injected (Fig. 7).
Note that the choice of a 2D problem was driven by numerical efficiency, as a 3D problem would significantly increase the model run time and prolong the planned sensitivity analysis.However, in light of the results shown in Section 3.2, we were quite satisfied with the accuracy of the 2D model.We also believe that a 2D model could be representative of the problem, considering that gas flow should be fundamentally radial.In future work, we aim to compare 2D and 3D problems to verify this assumption.
The simulation of gas flow injection is divided in 5 stages (Table 1): First, a preliminary stage, stage 0, used for the generation of the initial conditions.Then, stage 1, in which the drilling of the borehole is modelled.After that, the test interval is installed at stage 2. The intervals are initially saturated with water and then, a water/gas exchange is modelled at stage 3.The process culminates with a gas injection stage (stage 4), and finally, a recovery stage (stage 5).
In the PGZ1003 test, the injection stage was divided into two phases, each with an approximate rate of 500 mLn/min (1 mLn ≈ 2.2 ×10 -8 kg/ s), following the in-situ test protocol of Fig. 8.The decision to adopt this two-phase approach was driven by the pressure in the interval nearing the packers' pressure (approximately 14 MPa) after 90 min of injection.Consequently, the injection process was paused, the packer pressure was increased to 15 MPa, and the injection was resumed for an additional 35 min.
After gas flow finally ceases, the behavior of the model during the recovery phase is also observed.Interval #4, where gas injection occurs, was initially saturated with water.However, before beginning the gas injection test (stage 4), approximately 334 g of water were removed from the interval by applying a slight gas overpressure (stage 3).
In the model, the geometrical volume of the interval is 0.001115 m 3 , although the theoretical total volume of the interval is 0.001187 m 3 , accounting for the volume of the hydraulic lines leading to the control panel (0.000213 m 3 ), the volume of the technological void (0.000204 m 3 ), the volume of the filter (0.000111 m 3 ) and the volume of the existing annulus between the filter and the borehole wall (0.000659 m 3 ).However, due to convergence, the volume of the annulus could reduce to zero.Therefore, depending on the convergence of the borehole wall, the theoretical total volume of the interval may vary between 0.000528 m 3 and 0.001187 m 3 .Consequently, before starting the gas injection test, the degree of gas saturation in the interval may vary between 28% (without any convergence of the borehole wall) and 63% (full convergence of the borehole wall).
However, during the saturation of the intervals (stage 2), 0.001947 m 3 of water were injected to saturate the interval, resulting in a difference between calculated and actual volume injected of 0.000760 m 3 .To account for this discrepancy, a virtual porosity value is used which may vary between 1.06 and 1.75.Note that these virtual values of porosity greater than one have no physical meaning, but are necessary in this case to account for the total volume of water in the interval (e.g.porosity=2 means twice the geometrical volume in the model).For the base case, a virtual porosity value of 1.25 is considered for the interval.
Regarding boundary conditions, an anisotropic compressive stress state is considered, with a stress of 12.4 MPa applied horizontally and 12.7 MPa applied vertically at far boundaries.Normal displacements are not allowed on the left and lower boundaries.Before the installation of the interval, the displacements on the borehole's wall are not restricted.The effect of gravity is considered negligible for this problem.Moreover, a water pressure of 4.7 MPa and a gas pressure of 0.1 MPa are applied at the far boundaries, while left and lower boundaries are impermeable.
Regarding the initial and boundary conditions in the interval, which is installed in stage 2, initial porosity is considered to be equal to 1.25 and all initial pressures and stresses are considered atmospheric (0.1 MPa).Normal displacements are not allowed on the internal surface of the interval.During stage 2 the interval is saturated with water at a pressure of approximately 4.7 MPa.Then, during stage 3, a slight gas overpressure of around 4.8 MPa is applied in the interval.
As for initial conditions in the COx, the initial liquid pressure is 4.7 MPa, the initial gas pressure 0.1 MPa, and the initial stress state considered is anisotropic: 12.4 MPa in the x-axis, 12.7 MPa in the y-axis and 16.1 MPa in the z-axis.Temperature is considered constant at 20º C  The average hydraulic conductivity of the COx is approximately 4 × 10 -13 m/s, resulting in an average initial intrinsic permeability k 0 of approximately 4 × 10 -20 m 2 parallel and 1.33 × 10 -20 m 2 perpendicular to the horizontal bedding.The change in intrinsic permeability of the matrix with porosity has been modelled using the exponential law (Eq.4), with b k = 42.Initial porosity of the COx is considered to be 0.18 for the homogeneous base case.
In addition, several stochastic heterogeneous cases have been considered, each with a mean porosity of 0.18 and a variance of 0.00030.This results in an initial porosity range, excluding extreme values, approximately between 0.14 and 0.22.Fig. 9 illustrates three different stochastic distributions of porosity as an example.
In the heterogeneous cases, a small borehole damaged zone (BDZ) with a radius of 0.053 m is considered, as the porosity (and thus permeability) around the excavation should be higher in this zone due to the disturbance on the borehole's wall caused by drilling.This BDZ has a mean porosity of 0.22 and a small variance of 0.00005, yielding an initial porosity range approximately between 0.21 and 0.23.Note that during the creation of the heterogeneous porosity distribution, we aim to prevent the elements in the BDZ from starting with a low porosity value.
On the other hand, Andra has conducted high-resolution seismic surveys to better understand the distribution of the porous and permeable bodies for COx, 45 but the resolution of the 3D seismic grid is 10 m by 10 m horizontal and 6 m vertical. 46,47This resolution is not applicable to our model, which is a square of dimensions 10 m × 10 m.
In this article, the range values in the vertical direction (a y = 0.6 m) have been chosen to be bounded by sample size and element size, and the large values in x and z-directions (a x = a z = 30) have been chosen to force a marked horizontal anisotropy (Figs. 1 and 2).Finally, we have considered C 0 = 0 for simplicity.Hence, different stochastic anisotropic Fig. 9. Three distinct stochastic anisotropic initial distributions of porosity.Fig. 10.Variation in intrinsic permeability induced by porosity (left) and strain (right).Note that the left figure represents the initial intrinsic permeability k 0 values corresponding to the initial porosity distribution, according to the exponential law (Eq.4), whereas the right figure displays the variation in intrinsic permeability induced by strain (Eq.5), for an element with k 0 = 4 × 10 -20 m 2 .

Table 2
Hydraulic parameters for the COx in the base case of the 2D gas injection model.
Water retention curve (van Genuchten;Eq.6): Intrinsic permeability (Darcy's law;Eq.1): Gas phase relative permeability (Eq.7):Note that A is a gas enhancement factor attributed to changes in clay material structure.
initial distributions of porosity have been created (Fig. 9).The heterogeneous porosity field results in an initial intrinsic permeability of the matrix ranging from 8 ×10 − 21 m 2 and 2 × 10 − 19 m 2 , as depicted in Fig. 10 (left).This is calculated using the exponential law (Eq.4).In addition, a cubic law (Eq.5) has been employed for the embedded fracture model, as shown in Fig. 10 (right).The parameters for this model are calibrated and presented in Table 2.It is important to note that the porosity remains relatively constant throughout the experiment, due to the relatively high material stiffness and almost constant volume test conditions.Therefore, the exponential law primarily influences the initial intrinsic permeability, while the cubic law enhances the permeability as a function of the strain locally developed.Table 3.
The COx material behavior is characterized by anisotropic stiffness, with Young's modulus E of 6000 MPa parallel and 4000 MPa perpendicular to the horizontal bedding.Similarly, Poisson ratio υ is 0.3 parallel and 0.2 perpendicular to the horizontal bedding.The shear modulus is 1700 MPa, the uniaxial compression strength is 21 MPa, and the tensile strength is 3 MPa.The density of solid grains is considered to be 2770 kg/m 3 .In cases where a BDZ is considered, the parameters of the BDZ are identical to those in the COx, except for the initial porosity, as explained earlier.Regarding the interval, its density is considered to be 1631 kg/m 3 .Other interval parameters are listed in Table 4.

Results and sensitivity analyses of the short-term gas flow injection experiment
The evolution of gas injection pressure in the homogeneous base case model and in different heterogeneous models is shown in Fig. 11, and compared with experimental measurements.It can be observed that there is some variability in the gas pressure across the different cases, reflecting the variability that we can encounter in the field due to the heterogeneity of geomaterials.However, these differences are not significant and it appears that a homogeneous model can adequately represent the experiment.
However, modelling the gas injection test is complex and requires sensitivity analyses due to the uncertainty of some parameters.For instance, the degree of gas saturation S g in the interval before the injection can vary between 0.63 and 0.28, as commented above.Fig. 12 presents a sensitivity analysis to this parameter for S g = 0.4 (base case) and S g = 0.5.The results suggest that, although both curves represent well the behavior of the model in the recovery phase, a lower degree of gas saturation S g does not accurately represent the two pore pressure peaks caused by the two different injection phases during the gas injection stage (see Fig. 8).Therefore, we can conclude that the degree of gas saturation before injection must be around 0.4.
Another uncertainty, previously mentioned, is the volume of the injection chamber.Fig. 13 shows the pore pressure evolution of several homogeneous models with different interval volumes, ranging from the theoretical volume (1187 mL), corresponding to a virtual porosity of the interval of 1.06, to the actual injected volume (1947 mL), corresponding to a virtual porosity of 1.75.Intermediate cases with porosities equal to 1.2, 1.25 (base case) and 1.3 are also presented, corresponding to interval volumes of 1340 mL, 1400 mL (base case) and 1450 mL, respectively.It can be observed that results are similar in the recovery phase, like in the previous sensitivity analysis.However, the pore peak

Table 4
Parameters for the interval in the base case of the 2D gas injection model.pressures obtained are different (especially the first one).It can be concluded that considering the injected volume underestimates the pore peak pressures, while considering the theoretical volume slightly overestimates the pore peak pressures.Therefore, according to the numerical results, the volume of the interval must be in the range of 1400 mL, higher than the theoretical one but lower than the initially injected volume.
The influence of a BDZ with different hydraulic properties has also been studied (Fig. 14).To do that, a BDZ of radius = 0.053 m has been considered and compared to the base case (in which BDZ is not considered and, hence, porosity is 0.18 and intrinsic permeability is 4 × 10 − 20 ).This BDZ has a porosity of 0.22, which corresponds to an intrinsic permeability of 2.1 × 10 − 19 , approximately one order of magnitude larger than the intrinsic permeability in the COx.However, results indicate that considering a BDZ with the previously mentioned characteristics does not affect significantly the results when compared to a case without BDZ (Fig. 14).
It is important to note, as highlighted in Section 3.1, that all heterogeneous cases consider a BDZ with these specific characteristics.Therefore, we aimed to assess the impact of this particular BDZ in a homogeneous case to ensure that the disparities observed in the heterogeneous cases are not predominantly attributed to this BDZ.Nevertheless, it is important to clarify that this study did not examine the influence of a BDZ with different properties, such as those found in Armand et al. 48nother sensitivity analysis has been conducted to analyze the effect of plasticity and stiffness anisotropy (Fig. 15).To do that, the elastoplastic base case has been compared to an elastic case, with an additional elastic case considering material isotropy (with equivalent elastic parameters E = 5240MPa and υ = 0.3).According to the results, mechanical plasticity does not play an important role in modelling this particular short-term experiment as very few or none of the elements in the model reached their maximum strength.However, it is important to note that this may not hold true for other experiments where plasticity could potentially play a significant role, such as those simulated in Section 4. On the contrary, it is important to consider stiffness anisotropy, since the results are significantly different compared to the isotropic case (Fig. 15).
Lastly, Fig. 16 illustrates a sensitivity analysis to the embedded fractures spacing s value.Note that, while we have not experimentally measured this parameter, it does hold physical significance. 33,35ig. 17 depicts the shape of the different curves using different embedded fracture spacing (Eq.5): 5 × 10 -4 , 2 × 10 -4 , 1 × 10 -4 , 5 × 10 - 5 and 2 × 10 -5 m.As shown in Fig. 17, lower values of the spacing s reduce the slope of the change in intrinsic permeability with strain and consequently result in higher pore peak pressures.It has been found that this parameter has a significant influence on the results, being s = 2 × 10 -4 m (base case) the calibrated value for the modeled experiment.

Long-term gas flow injection modelling analysis
In Section 3, a real gas injection experiment has been modelled and calibrated, and sensitivity analyses have been performed on relevant parameters.However, the ultimate goal in this field of research is the long-term underground disposal of radioactive waste. 16,49Therefore, long-term analyses of gas flow injection have been conducted.In these    analyses, the recovery stage (stage 5) has been eliminated and a constant gas injection rate of 500 mLn/min is maintained instead.
In Fig. 18, different long-term analyses are represented.First, the effect of heterogeneity has been assessed.A homogeneous (base case) model, identical from the base case model in Section 3 but now maintaining the gas injection, has been compared to different stochastic heterogeneous cases also with the very same parameters, but heterogeneous (see Fig. 9).A heterogeneous elastic case has also been introduced.In view of the results (Fig. 18), a number of conclusions can be drawn.First, the maximum pore pressure obtained (the so-called breakthrough pressure, around 18.5 MPa), is higher than that obtained in the experiment (approximately 4 MPa higher).This suggests that the breakthrough pressure may not have been reached in the experiment.Moreover, the break-through pressure of the homogeneous case is slightly higher.
On the other hand, the results indicate that the consideration of heterogeneity or plasticity may not significantly alter the long-term outcomes.However, minor differences can be observed in the build-up of pore pressure across different stochastic heterogeneous models.
Another sensitivity analysis has been conducted on the embedded fractures spacing value, although this time focusing on the gas pathway (Fig. 19).It can be observed that different values of the embedded fractures spacing (see Fig. 17) result in different gas pathways.A relevant finding is that no hydraulic fracture develops in this model using a spacing value greater than approximately 5 × 10 -5 m.For greater values of spacing, the gas pressure develops without creating a preferential path.Moreover, preferential paths seem to be different for different values of spacing.In addition, plasticity is also necessary in order to obtain a preferential path (see the elastic case in Fig. 19).Therefore, a preferential path cannot be obtained in elastic models or in models with a spacing greater than a certain value (near 5 ×10 -5 m in this case).
Fig. 20 illustrates the long-term pore pressure evolution for homogeneous elasto-plastic cases, considering different embedded fractures spacing s values: 2 × 10 -4 (base case), 5 × 10 -5 , and 2 × 10 -5 m.As previously observed in Fig. 16, lower spacing values lead to higher break-through pressures due to the slower rate of change in intrinsic permeability with strain (see Fig. 17).
Moreover, an elastic case corresponding to s = 5 × 10 -5 m is also included in Fig. 20 for comparison.When comparing this purely elastic Fig. 18.Long-term pore pressure evolution within interval #4 (PGZ1003).Comparison between the homogeneous base case and various stochastic heterogeneous cases, including a heterogeneous elastic case.Experimental results corresponding to the 'short-term' gas injection in-situ experiment (discussed in Section 3) are included for the purpose of contextual reference only.case with the corresponding elasto-plastic case, significant differences emerge in terms of pore pressure evolution.Notably, the pressure breakthrough is nearly 2 MPa higher in the elastic case.This suggests that considering plasticity becomes relevant for spacing values lower than 5 × 10 -5 m, i.e. when a hydraulic fracture is created.However, its relevance may be lower when a hydraulic fracture is not developed.These differences between elasto-plastic and elastic cases become even more pronounced when considering a greater time span in the analysis, as shown in the lower graph of Fig. 20.
Finally, a few heterogeneous elasto-plastic models with a spacing s = 5 × 10 -5 m have been analysed (Fig. 21).Unlike the models in Fig. 18, in this case the spacing value allows the development of a hydraulic fracture.It can be observed that different initial stochastic distributions of porosity result in significantly different pressure evolutions when a hydraulic fracture is developed.The path of the hydraulic fracture developed also varies (Fig. 22) and tends to be horizontal due to the effect of horizontal bedding as reproduced by the anisotropic range used (Fig. 9).Therefore, considering heterogeneous models is relevant when a hydraulic fracture is developed, which in turn depends on the spacing value s.

Model validation with a similar experimental test
An experiment similar to that performed in borehole PGZ1003, discussed above, was conducted in a borehole named PGZ1002 (drilled in February 2020), oriented parallel to the minor horizontal stress (i.e.perpendicular to borehole PGZ1003).Most parameters of this experiment are identical to the previous one.Model geometry, mesh, materials, and boundary conditions are the same as in Fig. 7, although in this case, the applied stresses are 16.1 MPa in the horizontal direction and 12.7 MPa in the vertical direction.A slightly higher water pressure of 5 MPa was applied on the far boundaries for this test.
The stages of this test are the same as in the previous one (Table 1), although in this case the injection protocol has only one injection phase.The gas injection was performed in interval #4, which was initially saturated with water, although 472 g of water were removed from the interval by imposing a slight gas overpressure.In this test, the theoretical total volume of the interval is 0.001180 m 3 , accounting for the volume of the hydraulic lines leading to the control panel Fig. 20.Long-term pore pressure evolution in the interval #4 (PGZ1003).Homogeneous elasto-plastic cases considering different embedded fractures spacing s values: 2 × 10 -4 (base case), 5 × 10 -5 and 2 × 10 -5 m.An elastic case corresponding to s = 5 × 10 -5 m has also been added.Lower graph shows a greater time span for the same models.Experimental results corresponding to the 'short-term' gas injection in-situ experiment (discussed in Section 3) are included for the purpose of contextual reference only.

21.
Long-term pore pressure evolution in the interval #4 (PGZ1003).Heterogeneous elasto-plastic cases considering an embedded fracture spacing value s = 5 × 10 -5 m. (0.000206 m 3 ).The degree of gas saturation in the interval before gas injection may vary between 40% (without any convergence of the borehole wall) and 91% (full convergence of the borehole wall) for this case.
However, the volume of water injected to saturate the interval during the experiment was 0.001421 m 3 .Thus, there exists a difference between calculated and actual volume injected of 0.000241 m 3 .Hence, to match the total volume of water at saturated conditions, the value of porosity of the interval may vary between 1.06 and 1.27.The porosity value of the interval is considered to be 1.25; the same as in the test previously modelled (PGZ1003).
Fig. 23 shows the results regarding pore pressure evolution for the PGZ1002 test.With the gas saturation degree S g = 0.4 (the same as in PGZ1003 base case), it can be observed that the results compare well with the experimental data.The results of a case considering heterogeneity have also been added to Fig. 23.It can be observed that heterogeneity does not significantly change the results because no hydraulic fracture has been generated.
Finally, different long-term models with hydraulic fracture have been run to analyze the effect of heterogeneity (Fig. 24).The value for the embedded fractures spacing has been changed to s = 2 × 10 -5 m, so that a hydraulic fracture develops in the model.It can be noted that different stochastic heterogeneous distributions have a significant influence on the results, with differences higher than 2 MPa in the computed break-through pressure.

Conclusions
This article evaluates the influence of considering heterogeneous material properties on the results of a gas injection experiment performed in COx within a two-phase flow coupled hydro-mechanical (HM) framework.The approach used in the study offers a more realistic representation of subsurface conditions by accounting for the spatial variability of material properties using a variogram.The models have been executed in the CODE_BRIGHT program, a valuable tool for simulating coupled hydro-mechanical problems in heterogeneous geomaterials.
An experimental test (PGZ1003) performed in COx has been modeled and sensitivity analyses on different parameters have been performed, particularly those that are critical or uncertain.The results demonstrate that the experiment can be modeled successfully, as the experimental measurements have been accurately replicated.However, modelling the gas injection test is complex and requires sensitivity analyses due to the uncertainty of some parameters; in particular, the degree of gas saturation in the interval before injection, the volume of the interval, and the embedded fractures spacing require calibration.
Furthermore, long-term gas injection tests have been modeled.The results indicate that heterogeneity significantly impacts the results when a hydraulic fracture is generated.To generate a hydraulic fracture, it is necessary to consider plasticity and a small value of embedded fractures spacing.Note that, during the short-term gas injection experiment, a hydraulic fracture has not been generated.The reason for this is twofold: firstly, the volume of gas injected was insufficient to reach the breakthrough pressure; secondly, the calibrated spacing of the embedded fractures is not low enough.As a result, both the heterogeneous and homogeneous results appear similar.Finally, a validation exercise has been conducted in a different test (PGZ1002) using the calibration from the previous test performed nearby.The results align well with experimental data, validating the calibration of the model.
The study underscores that considering heterogeneity has, in certain conditions, a significant impact on the gas flow path and the breakthrough pressure.This article provides a framework for future research on the impact of heterogeneity on subsurface processes and the development of more precise and reliable simulation models.The study emphasizes the importance of considering the heterogeneity of material properties in the design and analysis of subsurface engineering projects.

Declaration of Competing Interest
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.

α non-associativity parameter that ranges from 0 to 1 α
, γ, β The angle with the positive direction of the coordinate axis (counterclockwise is positive) β parameter in Drucker-Prager yield function γstrain in the cubic law for calibration ε e , ε i elastic and inelastic strain tensor ε e v , ε e d elastic volumetric and deviatoric strain λ 1 range on the y axis divided by range on the x axis λ 2 range on the z axis divided by range on the x axis λ c σ ′ x , σ ′ y , σ ′ z effective stress in the x, y and z axis direction τ T vector of normally distributed numbers, transpose of the vector T Temperature A. Rodriguez-Dono et al.

Fig. 2 .Fig. 3 .
Fig. 2. Parameters required to define the geometric anisotropy of a covariance function in 3D.

Fig. 5 .
Fig. 5. Series of parallel embedded fractures with aperture b and spacing s.(adapted from Olivella et al.33 ).

Fig. 6 .
Fig. 6.Schematic view of the interval installed in the borehole.

Fig. 12 .
Fig. 12. Pore pressure evolution in the interval #4 (PGZ1003) compared with experimental results.Homogeneous cases with different initial gas saturation degree in the interval S g : 0.4 (base case) and 0.5.

Fig. 13 .
Fig. 13.Pore pressure evolution in the interval #4 (PGZ1003) compared with experimental results.Homogeneous cases with different values of the interval volume: from theoretical (1187 mL) to actual injected volume (1947 mL).

Fig. 24 .
Fig. 24.Long-term pore pressure evolution in the interval #4 (PGZ1002).Heterogeneous elasto-plastic cases considering an embedded fracture spacing value s = 2 × 10 -5 m.Experimental results corresponding to the 'short-term' gas injection in-situ experiment (discussed in Section 3) are included for the purpose of contextual reference only. 28

Table 1
Gas injection test timeline for interval #4 of borehole PGZ1003.

Table 3
Mechanical properties for the COx in the base case of the 2D gas injection model.