Effect of Compaction Banding on the Hydraulic Properties of Porous Rock - Part II: Constitutive Description and Numerical Simulations

In this paper the performance of a constitutive model for the description of the hydro mechanical behaviour of soft rock is evaluated with respect to the experimentally observed behaviour of Maastricht Calcarenite under different stress states that is presented in the companion paper. The mechanical model is elasto-plastic and consists of an associated yield surface, internal variables for the description of the hardening and softening behaviour and a non-local extension for the simulation of strain localization in form of shear bands and compaction bands. The model is implemented in the software ABAQUS and the laboratory results from the tests under dry condition with Maastricht Calcarenite are used for the calibration. The good agreement of the numerical results with the laboratory results is shown and the suitability of the model is discussed. To describe the effect of compaction bands on the permeability of soft rocks a simple analytical model based on the Kozeny–Carman equation is proposed and calibrated with the experimental results from drained tests under different stress states for Maastricht Calcarenite rock material. As the results are in good accordance with the experimental results, the model is implemented in the software ABAQUS and the numerical results are presented and discussed. Finally the performance of the model is evaluated and possible improvements are suggested.


Introduction
Strain localization in porous rocks like the formation of compaction bands, has been the subject of several studies throughout the last decades, including experimental and constitutive descriptions. On the experimental side, they have been investigated in the field, for example in Aztec Sandstone in Nevada (Sternlof 2006), in Navajo Sandstone in Utah (Antonellini et al. 1994) and in the Orange area in France (Ballas et al. 2012). But the formation of compaction bands has also been studied in the laboratory as in the companion paper. On the constitutive side, different theoretical approaches are applied to describe the onset of compaction bands. Olsson (1999) was the first that suggested the bifurcation theory for the description of compaction bands. Issen and Rudnicki (2000) adapted the analysis of shear bands from Rudnicki and Rice (1975) to compaction bands and demonstrated that compaction bands occur in porous materials that are described by a yield surface with a cap. Numerical analyses with the finite element method (FEM) (Das and Buscarnera 2014;Chemenda 2009) as well as the discrete element method (DEM) (Wu et al. 2020;Katsman et al. 2005;Dattola et al. 2014) were used to investigate the mechanism leading to the formation of compaction bands in porous rocks. The behaviour of bonded geo materials was studied via theoretical investigations (Tengattini et al. 2014;Nova et al. 2003). However, the calculations based on the assumption that the global behaviour and the local behaviour coincides although the laboratory results showed the opposite. To take the mechanical behavior on the micro scale into account non-simple-continuum models are necessary. The problem then does not show mesh-dependence in the simulations. In this study the application of this model is validated with experimental results of triaxial tests and oedometric tests on Maastricht Calcarenite, a soft rock with a weak cementation and a porosity of 52 % (Leuthold et al. 2020).
A second focus in this study is the constitutive description of the effect of compaction bands on the hydraulic properties of the material. In several laboratory studies with porous grain stones (Holcomb and Olsson 2003;Baxevanis et al. 2006;Fortin et al. 2005;Baud et al. 2012) a permeability decrease due to the formation of compaction bands was documented. The results from the drained triaxial and oedometric tests with Maastricht Calcarenite (Leuthold et al. 2020), demonstrate a decrease of the permeability of up to 3 orders of magnitude due the formation of compaction bands. Hence, one objective of this study was to establish a relation between the compaction process and the decrease of permeability in the rock. A common approach to describe the relation between void ratio and permeability is the Kozeny (1927); Carman (1997) equation. Ashari et al. (2016) as well as Nguyen and Einav (2009) proposed a continuum modeling approach where the Kozeny-Carman equation is combined with the breakage mechanics framework (Einav 2007 andTengattini et al. (2014)) to describe the effect of compaction on the hydraulic properties in porous rocks. In this paper, an analytical approach is presented where the interplay between the permeability and inelastic compaction is described with the Kozeny-Carman equation and the formation of compaction bands is considered.

Constitutive Description and Numerical Modelling of the Mechanical Behaviour
In the present section constitutive descriptions and the results of numerical simulations are presented for the mechanical response of the material. The experimental data is used to calibrate the parameters of a constitutive model reproducing the mechanical behaviour of the rock.
As an elasto-plastic model is used for the constitutive description of the material, a suitable formulation for the yield surface and the internal variables is needed. For the description of the yield surface a model by Gerolymatou (2017) is used. The yield surface is described by the following equation in the p-q stress space: where the multiplier M f is a function of the lode angle : The dependence of is expressed by the following function: (1) after Gudehus (1973), where c is defined: with M e as the inclination of the critical state line in extension and M c in compression. The function h(p) of the mean pressure p controls the shape of the yield surface: The parameters and are constant and determine the yield surface shape. The parameter p t is the tensile strength and p c the strength at isotropic compression. In Fig. 1 (Nova et al. 2003). In this model the scalar quantities p m , p s and p t Fig. 1 The experimental results of the triaxial tests under drained and under dry conditions and the calculated yield surfaces are introduced as internal variables to describe the loading history of the material. It is assumed that : where p s accounts for the effect of irreversible fabric modification and p m the effect of interparticle bonding. According to Nova et al. (2003) it is assumed that: with ̇p v as volumetric strain increment, ̇p q as deviatoric strain increment and the parameters s , m , s and m as material constants. The tensile strength is considered proportional to p m : where m is a dimensionless parameter. The constitutive equations are implemented in user defined subroutines, which can be used in the finite element (6) p c = p s + p m , software ABAQUS. The subroutines in ABAQUS are named UMAT (User MATerial) and are written in Fortran. For the calibration of the parameters the incremental driver code by Niemunis (2017) was used. The numerical calculations are implicit and only performed with one Gauss point. In Fig. 2 the experimental results of the triaxial tests under dry conditions and the results of the calculations with one Gauss point are presented. The numerical results are in good accordance with the experimental results, especially for high confining stresses. The hardening behaviour at the end of the tests and also the softening at the onset of plastic deformation can be described with the constitutive model. The values for parameters that were used in the simulations are listed in Table 2.
For the numerical simulation of strain localization as a boundary value problem with more elements, it is necessary to use a regularization, so that the results are not mesh dependent and the size of the deformation zone is realistic. To this end a non-local model was used. In this model the internal length l c is introduced. A bell shaped function was used:  The non-local average of the quantity u is defined then as a weighted integral of the local quantity as The internal variables p s and p m are averaged in the manner described above, resulting in a non-local expression for the strength. Figure 3 shows the averaging function for different values of the internal length.
Besides the weighted averaging of the internal variables, a viscosity is introduced in the material behaviour. For the sake of simplicity the Perzyna model (Perzyna 1966) was used here. The strain rates are decomposed as follows, while the evolution of the visco-plastic strain is defined as where g is the plastic multiplier, (f ) is the overstress function here assumed to be given by f ∕p c (t = 0) , f the current version of the yield function, is a normalized viscosity parameter given in sec and Furthermore, the visco-plastic strain rate evolves via a flow rule, with ̇ as the consistency parameter. The combination of the two equations (15) and (13)  This results in the elasto-visco-plastic modulus being expressed in the following equation: For the needs of the algorithm, the approximation is used.
In the next step plane stress tests with different side pressures were calculated using the software ABAQUS standard with an explicit time integration. In plane stress tests C3D8R elements are used and the CAE-model is only one element thick. The stress in the x1-and x2-direction is constant. The parameters that were used in the calculations are presented in Table 3.
To ensure that the viscosity has no influence on the numerical results, the biaxial test with 4 MPa has been calculated with different values for the viscosity . The reason for this choice is that no tests were available for rate dependence. The value of the viscosity was set to 100 s for all numerical calculations, because the results of the ratedependent calculations were similar to those of the rate-independent calculations and the calculation were more stable. In the present calculations the internal length is only used like a fitting parameter and has no physical meaning. However, from experimental observations the size of the deformation zone is known and a suitable value for the internal length would be the width of a discrete compaction band ( l c = 1 mm). As the element length needs to be 10 times greater than the value for the internal length an element length of 0.1 mm is necessary. This would lead to a high computation time so that it was chosen to set the internal length to 10 mm. The non-local model was implemented to alleviate mesh dependence in the numerical simulations. In the following, the influence of the element length on the results will be investigated. A biaxial test at a confining pressure of 3.25 MPa was calculated with an internal length of 10 mm and an element length of 1.4 mm and 0.7 mm. In Fig. 4 the distribution of void ratio is plotted for the same value of the axial strain for different values of the element length. The size of the shear band is independent of the used mesh and shows the same dimensions for both simulations. That means the non-local model is suited to alleviate mesh dependence and to determine the size of the deformation zone through the internal length.
The test duration was set to 10800 s, equal to the test duration in the laboratory and the time step was set to 0.15 s. To control the location of the onset of the localization, an initial defect was introduced in the samples by reducing p m slightly in the range of 5 mm in the middle of the sample. In plane stress tests all stresses in the horizontal direction, including the out of plane ones, were set equal to the confining pressure and maintained constant. In Fig. 3 the results of the simulations with 1 MPa, 2.5 MPa and 4 MPa confining pressure are shown.
The distribution of the void ratio is presented in Fig. 5. For high confining pressures at 4 MPa the sample fails by the formation of a compaction band and for confining pressure at 1 MPa, the sample fails by the formation of a shear band. The mechanical behaviour of the rock is similar to the behaviour observed in the experiments. However, instead of the formation of discrete compaction bands, the compaction process is homogeneous and when the process is finished the void ratio in the whole sample is uniform.
In Fig. 6 the axial stress is plotted against the axial strain for the biaxial tests with different confining pressures and the corresponding experiments.
Compared to the results from the calculations with the material point, the mechanical behaviour is more brittle for low confining stresses up to 3.25 MPa. The major difference between the calculations is the non-uniform strain field in the biaxial tests. However, the numerical results are in good agreement with the experimental results, at least for the compaction. The spontaneous unloadings in the curve from the experimental test with 4 MPa confining pressure are a result of the formation of discrete compaction bands and the specimen inhomogeneity and therefore absent in the simulation. For low confining stresses like 1 MPa the  Das and Buscarnera (2014) this parameter has an effect on the width of the so called 'localization zone' in the yield surface and hence higher values for the parameter could lead to more suitable results. Figure 7 shows the experimental and the numerical results for the circumferential strain evolution as a function of the axial strain at different confining pressures. The numerical results are qualitatively and quantitatively in good accordance with the experimental results at small and at high confining pressures, while some deviation is present in the intermediate confining stresses. In Fig. 8 the results of a simulation of a triaxial test with a confining pressure of 4 MPa are plotted. The spatial distribution of void ratio is plotted for different values for axial strain. It is clear that the sample exhibits the formation of a compaction band and with further axial deformation the compaction zone propagates to the ends of the sample. This is not in agreement to what was observed in the experiments. However, the numerical specimen was uniform, something that is not the case for the actual specimens. Of course, it is also possible that the qualitative difference may arise from the calibration. In order to investigate the influence of the value of the internal length, biaxial tests with an element length of 0.7 mm and an internal length of 5 mm have been calculated. In Fig. 9 the differential stress is plotted against the axial strain for different confining pressures. For high confining pressures the numerical calculation with an internal length of 5 mm aborted at smaller values for the axial strain. In the test at 4 MPa confining pressure the calculations aborted before the compaction was completed. This is in agreement with the conclusions of Dattola et al. (2015) that higher ratios of specimen length to internal length are more unstable. A   different effect can be observed in the results of the test with small confining pressures like 1 and 1.5 MPa. In these tests the results do not coincide for the different internal lengths and the softening behaviour is more pronounced compared to the results of the calculations with an internal length of 10 mm.

Constitutive Description and Numerical Modelling of the Hydraulic Behavior
In the present section the effect of compaction banding on the permeability in axial direction is investigated with an analytical approach and numerical simulations.

Permeability Evolution: Analytical Approach
In the analytical approach the permeability was calculated with the Kozeny (1927); Carman (1997) equation. In this equation the permeability is described in terms of the void ratio and the specific surface: where K is the permeability, C kc is the Kozeny-Carman factor, S is the specific surface and e is the void ratio. Experimental results by Baxevanis et al. (2006) and Papazoglou et al. (2017) showed that at high confining pressures the specimens fail by the formation of compaction bands, meaning the deformation behavior is not uniform. Therefore, the compaction process needs to be considered in the calculation of the permeability. If the lateral deformation of a sample is neglected or inhibited, the permeability can be described as a function of the axial strain. During the compaction process, the sample consists of compacted and uncompacted parts. As a result, the height of the sample can be divided in a height for the compacted rock and one for the uncompacted rock.
The height for the compacted part is h c and the height for the uncompacted part is h u . In Olsson and Holcomb (2000) the mass balance is used for the derivation of the following equation for the height of the compacted part: In this equation e c is the void ratio of the compacted rock, e 0 the initial void ratio of the uncompacted rock, h 0 the initial height and pl ax is the plastic axial strain. The Kozeny-Carman pl ax equation is used to derive the following equation for the permeability of the compacted rock: In this formulation S 0 is the initial specific surface, S c is the specific surface of the compacted rock mass, K 0 is the initial permeability and K c is the permeability of the compacted rock. The ratio of the two specific surfaces can be seen as a grain crushing factor, because the specific surface increases with increasing grain crushing. The factor can be determined by micro structural analysis or by the determination of the permeability of the compacted rock, assuming that the Kozeny-Carman equation holds. The compacted and the uncompacted area may be viewed as connected in series, as far as the permeability is concerned. This yields the following relationship for the permeability of the specimen as a function of the permeabilities of the compacted and the uncompacted material The influence of the radial strain due to the change in crosssection has not been included here, as it was not measured in the present study. Its consideration is however algebraically straightforward and this effect could be considerable in the presence of large deformations. The simple analytical model predicts an anisotropic permeability parallel and perpendicular to the compaction bands. An example is given in Fig. 10, where the predicted vertical and horizontal permeability are given, as well as the prediction according to Kozeny-Carman. The initial and final void ratios were 1.0 and 0.6 respectively and both radial strains and potential grain crushing were ignored. In the vertical direction the consideration of (25) Fig. 10 Example of predicted permeability evolution the non-uniformity in the sample results in a more abrupt decrease initially.
The comparison of the predicted to the experimental results is difficult in the present case, as the radial strains and the void ratio of the fully compacted material are not known. The void ratio was chosen to fit the final void ratio to the results of each experiment. Grain crushing was considered in the change of the specific surface as a fitting parameter, in the absence of measurements. The parameters for the model calculation are presented in Table 3. In Fig. 11 the results of the model and the experimental results are plotted for different confining pressures. The permeability in logarithmic scale is plotted against the axial strain. While the agreement with the experimental results is good, it is important to validate this simplistic model in the future against experimental data for which all parameters are known, to truly assess its suitability. The values for the void ratio and the relation of the specific surfaces S 2 0 S 2 c , the so called grain crushing factor, increase with increasing confining stress. This is reasonable, because the experimental results from triaxial tests under dry conditions show that high values for the confining stress lead to compaction in radial direction. Moreover, visual inspection of the specimen after drained triaxial tests at a confining pressure of 4 MPa and 5 MPa documented pronounced deformation in radial direction. Papazoglou et al. (2017) investigated the effect of compaction in triaxial stress conditions on the grain size distribution. He compared the grain size distribution of the intact material and the compacted material after a test at a confining pressure of 4 MPa and he showed that the proportion of fine grain increases due to compaction. It is however clear that from a qualitative point of view it is more compatible with the experimental results than the standard Kozeny-Carman formulation, which neglects the non-uniformity of the specimen. Figure 12 shows experimental results of the oedometer test under drained conditions and results of the calculations with the analytical model. The results of the calculation are in good accordance in the first part of the test, but at an axial strain greater than 20% the permeability is overestimated. Moreover, the curve shape of the experimental results is different compared to the results of the analytical model from this point, which also coincides with the onset Fig. 11 Results of the analytical model and experimental results from triaxial tests for different confining pressures: permeability in axial direction plotted in logarithmic scale against axial strain. a 2.5 MPa, b 3.25 MPa, c 4 MPa d 5 MPa 1 3 of hardening. Beyond this point the model underestimates the permeability decrease with increasing axial strain. In the thin sections grain crushing was documented, so that one possible reason for the difference of the results could be the increase of the specific surface. Another option is that the Kozeny-Carman equation is not suitable to describe the permeability evolution during compaction of porous rocks in oedometer tests. This is in agreement with experimental results by David et al. (1994) where the permeability was measured at different values for the porosity during isotropic compression tests with four different porous quartz sandstones. The evolution of the permeability as a function of void ratio is described well with a power law and the values for the exponent are between 4.6 and 25.4. This indicates that also for these results the Kozeny-Carman equation is not suitable to describe the hydraulic behavior. Another reason for this could be that the turtuosity is not considered in the Kozeny-Carman equation (Table 4).

Permeability Evolution: Numerical Simulation
To further investigate the influence of the non-uniformity resulting from strain localization on the permeability evolution, the Kozeny-Carman equation was used to describe the local permeability evolution in simulations of oedometric tests under drained conditions with the software ABAQUS. The parameters used in the simulations are listed in Table 3.
For the numerical simulations with ABAQUS it was necessary to calibrate the internal variables for the softening and for the hardening behaviour. The calibration was also performed with the incremental driver by Niemunis (2017). The parameters resulting from the parameter calibration are listed in Tables 1 and 2. In Fig. 13 the axial stress is plotted as a function of axial strain for the experimental results and the results of the parameter calibration. For the initial stress conditions of the oedoemeter test, the heating of the ring has to be considered. It was assumed that the sample was pre-stressed in the radial direction due to the cooling of the steel oedometric ring, which had an initial temperature of 200 o C. This results in an initial radial stress equal to 0.8 MPa. The results of an oedometer test simulation under drained conditions are presented in Fig. 14. On the left is the permeability plotted in logarithmic scale against the axial strain and on the right the axial stress is plotted against axial strain. The curve propagation of the axial stress as a function of axial strain for the experimental results is different to the one of the simulation. The main difference is that in the laboratory test the hardening starts at an axial strain of 17 % and in the simulation the inclination of the curve does not change during the plastic deformation. The permeability decrease due to the compaction of the sample is underestimated in the simulations, because grain crushing and its influence on the specific surface is not considered. To overcome these problems it is necessary to implement the change   Fig. 13 Axial stress versus axial strain for experimental results and the results from the parameter calibration at different confining pressures under drained conditions of the specific surface due to grain crushing in the numerical model. This could be done by describing the specific surface as a function of the stress state. As the formation of discrete compaction bands cannot be reproduced in the numerical calculations and the compaction begins in the middle of the sample, the pore pressure distribution in the sample did not change during the test.

Conclusions
In this study the formation of compaction bands and their effect on the hydraulic properties were studied by constitutive description and numerical simulations. For the constitutive description a non-local elasto-plastic model is used and validated with experimental results (Leuthold et al. 2020). The formulation of the yield surface and the internal variables are suitable for the description of the mechanical behaviour of soft porous rocks. The mechanical response was simulated both at the element level and as a boundary value problem, using a non-local formulation to alleviate mesh dependency. In the simulations the brittle ductile transition of the rock was in good accordance with the experimental results. However, in the simulations no discrete compaction bands were observed, but rather a continuous propagation compaction front. This may be attributed to the fact that the numerical specimen was uniform or may be a matter of parameter calibration. The characteristic length used for the simulations was found to have very little effect on the mechanical response in the domain of ductile failure, but affected the results when shear bands were observed. The difference between global and local behaviour should therefore be considered, when calibrating constitutive models.
For the evolution of the permeability a simple model using the Kozeny-Carman equation was developed, taking into account the non uniformity of the samples. It was found to fit the results significantly better than the standard Kozeny-Carman formulation, especially from the qualitative point of view. This seems to indicate that the sample non-uniformity is the reason for the deviation of the results from the response predicted by the Kozeny-Carman formulation, as observed by other authors already at small strains (Baxevanis et al. 2006). Under oedometric conditions it was found that the permeability is underestimated when grain crushing is not considered in the model.
Finite element simulations of drained oedometric tests showed that the permeability evolution is also underestimated in the simulations. As already mentioned, it is possible that this problem could be solved by linking the specific surface with the stress state. This indicates that the nonuniform strain field is not the sole cause of the deviation from the Kozeny-Carman law.
Funding Open Access funding enabled and organized by Projekt DEAL.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. Fig. 14 Results of the simulation and experimental results for the oedometer test under drained conditions: a permeability in the axial direction plotted in logarithmic scale against axial strain b axial stress plotted versus axial strain