{'fraction_non_alphanumeric': 0.024569186145708924, 'fraction_numerical': 0.0065176591025422285, 'mean_word_length': 5.007790077900779, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Continuum modeling of the effect of surface area growth due to crushing and damage on the permeability of granular rocks

This paper discusses a continuum approach to track the evolution of permeability in granular rocks by accounting for the combined effect of porosity changes, grain breakage and cement bond damage. To account for such a broad range of microscopic processes under general loading paths, the Breakage Mechanics theory is used and the computed mechanical response is linked with the Kozeny equation, i.e. a permeability model able to evaluate the reduction of the hydraulic conductivity resulting from the simultaneous loss of porosity and growth of surface area. In particular, the evolution of the internal variables of the model has been linked to idealized geometric schemes at particle scale, with the goal to distinguish the contribution of the fines generated by the disaggregation of the cement matrix from that of the broken fragments resulting from the crushing of the skeleton. Compression/flow experiments available in the literature for different granular rocks are used to validate the proposed methodology. The analyses illustrate that the drop of the permeability of damaged rocks would be severely underestimated without an accurate computation of the growth of surface area, as well as that the distributed fragmentation of skeleton particles tends to have stronger implications than the generation of cement fines. These findings, along with the satisfactory agreement between model predictions and experiments, stress the benefits of adopting microstructure-based constitutive laws for the analysis of coupled hydro-mechanical problems.


Introduction
Fluid flow in reservoir rocks is invariably coupled with various deformation mechanisms associated with grain breakage, compaction, strain localization and fracture. The interplay between such microstructural alterations of the solid matrix and pore-scale flow controls the overall permeability of rocks. Understanding such relations is pivotal for the success of any activity involving injection/extraction of fluids, such as oil and gas extraction and underground storage of by-products [1][2][3]. For granular porous media, the permeability is a function of the grain size distribution, as well as of how the degree of particle sorting evolves during deformation [4,5]. Such distributed fragmentation of the skeletal grains is often accompanied by the propagation of micro-cracks across the cement matrix bonding the grains, which in turn represents a further cause of surface area growth [6,7].
The variation of permeability in rocks has been widely studied from an experimental standpoint [8][9][10][11] and numerous empirical correlations were proposed [12,13]. Similarly, a wide range of numerical [14,15] and analytical studies [16,17] have been dedicated to this topic. The physical interaction between fluid flow and deformation processes; i.e. hydro-mechanical (HM) coupling, in rocks has also attracted considerable attention [18,19]; especially when flow-deformation analyses can be dramatically affected by the porosity-dependence of rock permeability. In this context, usual HM continuum formulations for granular solids tend to lack coordinated mechanical and hydraulic constitutive laws based on the same microstructural descriptors, and often resort to empirical permeability functions only coupled with porosity changes.
Several models, however, nowadays can take into account the evolution of the material microstructure during inelastic deformation, thus providing new avenues to account for the important effects of the changes in specific surface area that may result from micro-fracturing and grain crushing.
This paper aims to discuss a continuum modeling approach that quantifies changes in granular rocks permeability by tracking the evolution of microstructural descriptors linked with the simultaneous loss of porosity and growth of surface area. For this purpose, a continuum model originally developed based on the Breakage Mechanics theory [20][21][22] is used to track the progression of microscopic grain fragmentation and cement disaggregation through two dedicated internal variables, a Breakage index B and a Damage index D, linked with the grain size distribution and cement size distribution, respectively. Such micromechanical features, together with the ability to track changes in porosity and surface area, make the Breakage Mechanics theory suited for hydro-mechanical analyses at the constitutive level [23][24][25].
Former studies have indeed demonstrated its applicability to the study of permeability evolution in rocks by linking different classes of Breakage models to a Kozeny law describing the hydraulic conductivity of granular rocks [26,27]. In the following, the latter idea is further explored for the case of unconsolidated granular rocks (Glass Beads and Ottawa sand) and weakly cemented granular rocks (Otter Sherwood sandstone), with the goal to evaluate the performance of the Breakage Mechanics model against the available experimental data.

Permeability models for granular materials
Numerous empirical formulae, capillary models, statistical models and hydraulic radius theories have been proposed to estimate the permeability of granular rocks. Capillary tube models are among the simplest and most widely used models from which many other existing permeability laws are formulated, in that they provide a mechanistic explanation for the Darcy's law controlling the macroscopic flow of a viscous fluid through a porous medium. In a one-dimensional context such laws can be written as follows: where Q is the flowrate of the formation, κ is the intrinsic permeability of the medium, A is the crosssectional area available for flow, µ is the viscosity of the fluid and ∂P ∂x is the pressure gradient. The negative sign in Equation (1) indicates that the fluid flows from high pressure to low pressure. The capillary tubes model is formulated from the Hagen-Poisseuille's law which governs the steady flow in a single straight circular capillary tube of diameter d [28]: where ρ is the fluid density, g is the gravitational acceleration and ϕ is the pressure head. The analogy between Equation (1) and Equation (2) is apparent considering that ϕ = P/ρg for static pressure, while the medium permeability κ equals d 2 /32. The capillary tubes model considers a bundle of such parallel tubes, all of the same diameter d similar to Figure 1 and gives the specific dischargeq through the porous block for the case of N tubes per unit area of cross-section of the model as: Since the porosity n of this model is defined as n = N (πd 2 /4)L/L, Equation (3) can be rewritten as: which is again analogous to Darcy's law.
These principles are used in the Kozeny model [29], in that the porous medium is treated as a bundle of capillary tubes of equal length. The cross-section of the tubes is not necessarily circular in this model, but similarly to the original capillary tubes model only the fluid velocity component oriented along the tube direction is taken into account in the model derivation. By expressingq with respect to a unit volume of solid rather than a unit volume of the porous medium, the Kozeny model derives: where c 0 is the Kozeny's constant (e.g., c 0 = 0.5 for circular channels) and S is the specific surface area (i.e. the total surface area of solid material per unit of volume) which is inversely proportional to the hydraulic radius R of the porous medium (i.e. the ratio between the cross-sectional area and the wetted perimeter) as S = n/R. The expression of κ in Equation (5) is referred to as Kozeny's equation and has been used widely by many researchers to evaluate the permeability of granular media, and it will be used hereafter to assess the evolution of granular rock permeability in presence of simultaneous changes of porosity and surface area. .

Breakage Mechanics models
Modeling the deformation of geomaterials in a continuum mechanics framework implies challenges associated with the representation of the macroscopic effects of local inter-particle interactions. Within this context, crushing phenomena are characterized by the further challenge of simulating the implication of an evolving particle gradation.
The Breakage Mechanics framework [20,21] is an example of continuum theory for crushable granular materials in which the particle-continuum duality is captured by combining thermodynamics principles with micromechanical arguments. Specifically, the theory is capable of describing the evolution of the grain size distribution (gsd) by means of a breakage internal variable associated with the release of elastic strain energy due to local grain fracture. This feature creates an explicit link between micro-scale mechanisms and macro-scale mechanical response, thus offering the significant advantage of tracking microstructural attributes in addition to usual global measures of deformation.
This framework has been recently augmented by Tengattini et al. [22] to cope with cemented granular systems. Specifically, an internal variable reflecting the development of micro-cracks in the cement phase was added to the formulation, with the goal to mimic the implications of the progressive rupture of the bonds bridging the particles ( Figure 2). This has led to an enhanced Breakage-Damage formalism particularly convenient to track surface area changes in lightly cemented granular materials, where the damage of the cement interacts with the fragmentation of the skeletal particles.
Hereafter the components of the Breakage Mechanics formulations essential for this work are briefly discussed, while a more detailed description can be found in Einav [20] and Tengattini et al. [22].
According to the Breakage Mechanics theory, gsd can be used to define the current state of a crushable granular material and it can be tracked by means of a scalar internal variable named Breakage B, defined as: B = B t /B p . B t and B p can be deduced by using the initial and ultimate cumulative gsd functions F 0 (d) and F u (d) similar to Figure 3a where d represents the grain size. While the initial gsd is typically measured to define the initial state of the material, the ultimate gsd is often assumed to be represented by power-law distribution of fractal nature, similar to Equation (6) with a fractal dimension α [30,31].
whre d m and d M are the minimum and maximum grain size at the ultimate gsd and the fractal dimension α is often chosen in the range 2.5-2.9. The current cumulative gsd can be expressed as a function of B as follows: Rock grain Cement Broken grain Fines Figure 2: Broken cement between grains and the generated pseudo grains (fines).
For cemented granular materials, micro-cracking of the cement phase is associated with a decrease in the grain size (log scale) cement size (log scale) size of the mechanically active inter-granular bond, which results into cement bridges of smaller crosssectional area. Using this assumption, it is possible to define an effective cement size (x c ) distribution, which describes the evolution of the damage in the cement phase. Similar to breakage B, another scalar variable, i.e. damage D, is defined in the Breakage-Damage theory as: D = D t /D p . Again D t and D p can be determined by using the initial and ultimate cumulative effective cement size distribution Figure 3b. The initial effective cement size distribution can be calculated from X-Ray tomographic imaging of the material [32] while the ultimate effective cement size distribution can be hypothesized as a constant function implying the theoretical final state is a completely uncemented granular system; i.e., D = 1.
The current cumulative effective cement size distribution can be defined as The evolution of such two internal variables is linked with the storage and dissipation of energy, which implies the repartition of the work input done at the boundaries of an RVE,W , into the sum of the rate of stored energy (which under isothermal conditions can be quantified through the Helmholtz Free Energy,Ψ) and the rate of energy dissipation,Φ;W =Ψ +Φ [22]. The Helmholtz Free Energy in a system with multiple interacting solid phases can be expressed as the volume-average of strain energy contributions associated with the cement and granular phases, as follows: where ϕ g and ϕ c are the the volume fractions of grain and cement (ϕ g + ϕ c = 1). Ψ g and Ψ c are the average stored energy within the grain and cement volume fractions, respectively, which can be written as functions of the elastic strain energy in an unbroken/undamaged reference system (Ψ g r for the grain and Ψ c r for the cement phase): where ϑ g is a grading index emerging from a statistical homogenization procedure based on initial and ultimate grain size distributions [20]. The dependence of the Ψ function on breakage and damage indices emerges from statistical homogenization techniques, according to which the change in elastic properties is linked to the release of strain energy associated with permanent changes of the microstructural attributes.
The reference energy (for both grain and cement) can be subdivided into two distinct parts: the volumetric part Ψ v,g r , Ψ v,c r and the shear part Ψ s,g r , Ψ s,c r ; i.e. Ψ g r = Ψ v,g r +Ψ s,g r and Ψ c r = Ψ v,c r +Ψ s,c r , and each of these reference elastic potentials can be chosen to be either linear or nonlinear (pressure-dependent) functions.
Different definitions of these functions as well as the specific forms used in this work are presented in the Appendix.
Starting from the Helmholtz Free Energy, the Breakage Energy (E B = − ∂Ψ ∂B ) and the Damage Energy (E D = − ∂Ψ ∂D ) can be computed, which represent the strain energy released by the material during the breakage of the granular skeleton or the fragmentation of the cement bonds, respectively. Finally, by employing a coupled dissipation function based on the simultaneous presence of contributions due to breakage, damage and cohesive-frictional processes [21,22], it is possible to derive the following yielding condition: where M is the critical state parameter which is a function of the friction angle φ; i.e.
is the mean stress and q = ∂Ψ ∂ε e s is the differential stress. E BC and E DC are the critical breakage and damage energies; i.e. the energy thresholds at which the respective phases initiate the inelastic processes, in other words, material parameters to be calibrated.

Connection between the permeability model and the Breakage Mechanics models
As it can be seen in Equation (5), the permeability is a function of porosity, as well as of the surface area of the solid. Therefore for geomaterials exhibiting extensive micro-fracturing due to grain crushing and cement damage, the permeability tends to evolve not only because of porosity changes, but also because of alterations of their initial surface area [33]. While most constitutive models for granular materials do not incorporate the effects of an evolving microstructure into permeability predictions [34], breakage models can circumvent this problem by linking the evolution of the surface area to their internal variables B and D [35].
This logic can be illustrated in case of a simple assembly of spherical particles, for which the specific surface area S can be written in terms of the harmonic mean grain size D H . A specific grain fraction of the assembly, the surface area scales inversely with the hydraulic diameter, such as: S = 6/D H [28].
Hence, for a polydisperse material characterized by a probability density distribution for the grain size, the hydraulic diameter can be computed as the harmonic mean across all the fractions in the assembly, as follows: If a fractal distribution is considered for f (d), then D H will be calculated as: In the case of a crushable granular material, the evolving f (d) becomes a function of both B and d similar to Equation (7): By plugging this expression in Equation (12), the evolving harmonic mean grain size is derived as: where D H0 and D Hu are the harmonic mean grain size for the initial and ultimate probability density functions of grains (Equation (12)). The example above provides an expression of D H valid for a crushable assembly, and which can be readily computed as a function of the internal variable B. In the following, it will be shown that this logic can be extended to the case of cemented granular materials by using data about the accumulation of fines resulting from the disaggregation of a damaged cement matrix to derive an expression similar to Equation (14) in which breakage and damage contributions affect simultaneously the surface area growth.

Model simulations
In this section, the Breakage Mechanics model is calibrated and validated against the hydro-mechanical response of different granular solids. In other words, the model parameters are chosen to capture both the compression response and the evolution of the grain/cement size distributions under different loading paths and thereafter, the model is used to estimate the permeability reduction upon compression; i.e.
the ratio between the current permeability κ and the initial permeability κ 0 , by applying the Kozeny equation (Equation (15)) as follows: where n, n 0 are the current and initial porosity.

Ottawa sand
Ottawa sand is the first granular solid studied in this section. To start, the Breakage Mechanics model response under the oedometer stress path is calibrated against the experimental results reported for an initial porosity of 37% in [36]. The calibrated compressive response is plotted in (Figure 4a  Although a good agreement between data and model predictions is found, some discrepancies in the magnitude of permeability drop are noticed, with the model overestimating the permeability decrease compared to the experiment. Such mismatch tends to be larger at intermediate stress levels, and it can be related to the tendency of Equation (7) to overestimate the amount of breakage-induced grains during the first stages of post-yielding compression due to the assumption of using a linear scaling of the evolving grain size distribution. The mismatch, however, disappears as high compressive stresses are approached, i.e. when the computed gsd curve captures accurately the gsd measured at the end of the experiment.

Glass Beads
As a further example, here the variation of permeability in Glass Beads subjected to compression is discussed. Similar to Ottawa sand, the Breakage Mechanics model was calibrated to capture the mechanical response under the effect of radial loading paths. The latter were imposed by controlling a stress path parameter r defined as the ratio between radial and axial stress increments, r = ∆σ r /∆σ a (r = 1 representing hydrostatic paths). The experimental data for calibration and validation purposes are derived from previous studies reported by [38]. Figure     mental data and despite the limited dataset used for parameter calibration they provide a satisfactory estimate of the permeability at the end of the compression stage [38].

Cemented granular rocks
Compression processes may alter the microstructure of a granular solid in very different ways depending on whether or not a cementing agent is present. An example of this concept is illustrated in Figure   6, where the gsd of a cohesionless granular material before and after loading is compared with the initial and final gsd curves of a weakly cemented sandstone. It is readily apparent that, while for the granular material the increase in particle dispersion is accompanied by the preservation of unimodality, in the sandstone the compression process exacerbates the bimodality of the initial particle distribution. The second peak emerging at the lower bound of the distribution, indicates the generation of small grains that can be related to the disaggregation of the binding agent and the creation of fine particulate (hereafter more simply referred to as "fines"), (Figure 2). Since these fines participate to the evolution of the surface area and can have important implications on the hydraulic conductivity, (i.e. they clog the pores and increase the solid surface interacting with the fluid), they have to be taken into account in the calculation of the grain size distribution and consequently of D H . For the sake of simplicity, hereafter it will be hypothesized that all the fine particulates emerging from the disaggregation of the cement hinder the flow of water, thus further exacerbating the permeability loss. This hypothesis will provide a lower bound to the computed values of permeability, which will be compared against the scenarios in which the effect of the cement fines is entirely neglected.
Here, to estimate the accumulation of fines in relation with the damage of the cement matrix, the data a b Grain Size [7m] u (x f ). This relation implies that no fines are present when the cement matrix is undamaged, while the fine distribution achieves ultimate condition when the cement is locally fully damaged. As a result, the general expression for the grain size distribution considering both the skeletal particles and the fines generated by the disaggregating cement can be obtained as: where ϕ f is the current volume fraction of the fines calculated as ϕ f = Dϕ c . By analogy to Equation (14), the general expression for D H is given by: (17) which underpins an evolution of the material surface area as a function of both Breakage and Damage.
From Equation (17, it is readily apparent that by setting ϕ f = 0, (16) and (17) reduce to (7) and (14) respectively, thus recovering the relations valid for cohesionless granular systems as a particular case.

Otter Sherwood sandstone
Otter Sherwood sandstone is a weakly consolidated rock which has been characterized through hydromechanical laboratory tests by [39]. Since this material is weakly cemented, a cement volume fraction ϕ c = 10% has been used, in agreement with the range reported in the literature for rocks from the same formation [40]. The isotropic compression test is the first loading path used for model calibration, and it has allowed to obtain the model response shown in Figure 7a (calibrated parameters: K g =K c = 11900,  In addition, the permeability approximation using the D H for only grains and the combination of grains and fines are very close to each other as well as the values of permeability measured in the experiment and provide the upper and lower bound scenarios for permeability prediction depending on the amount of the fines participated in the permeability reduction processes. For this specific cemented granular rock, the prediction obtained by taking grain breakage into account and excluding the cement fines from the analysis is found to be the one closest to the experimental data, thus suggesting that the contribution of the fines to the overall permeability reduction tends to be much smaller than that associated with the crushing of the grain skeleton. From the same figure it is also readily apparent that the values of permeability computed by the models accounting for surface area growth are significantly more accurate than those based only on porosity changes, from which a considerable underestimation of the permeability losses emerges at all the stages of the simulation. While the hydrostatic test is used to define the compressive properties and the proportion between breakage and damage dissipation, information about the pressure dependence of yielding is used to define the properties associated with the shearing behavior, such as shear modulus, friction angle and cohesion (i.e., all properties that affect the size and shape of the yield envelope predicted by the model). For this purpose, the yield function (Equation (11)) is plotted in Figure 8    once again that the permeability reduction resulting from the predicted degree of grain breakage again proves to be the most accurate. By contrast, the estimates based only on porosity changes are found again to severely underestimate the permeability loss, while those accounting for the potential effect of cement fines tend to exacerbate the computed rate of permeability reduction.
The calibrated model can eventually be used to predict the mechanical response (Figure 10a), the grain size distribution at the end of the experiment (Figure 10b), the evolution of D H (Figure 10c) and of the permeability (Figure 10d) for the radial path loading of r= 0.6. By comparing the model prediction with the experimental data, a good agreement is observed for the mechanical response and the permeability reduction. A significant mismatch between data and model predictions can instead be observed in terms of final gsd for the test performed with r= 0.6. In the latter case, a much larger generation of fines was  Permeability approximation for Otter Sherwood sandstone under radial path loading of r= 0.6 [data re-plotted after [39]].

Conclusions
The analysis of permeability changes in granular rocks subjected to high compressive stresses has to account for significant microstructural alterations such as grain crushing, cement damage and pore collapse. To address this challenge, this paper uses the Breakage Mechanics framework and links an inelastic constitutive law for crushable granular solids to a permeability model for granular rocks. The proposed model tracks the changes in porosity and surface area by expressing them as a function of internal variables associated with the evolution of the grain size distribution and the cement size distribution of the rock. This choice enables the coupling of the deformation response with permeability evolution laws, and has been here exploited by calibrating the model parameters for two different granular materials (i.e., Glass Beads and Ottawa sand) and a weakly cemented granular rock (Otter Sherwood sandstone). The former set of analyses are aimed to study a simple proxy of unconsolidated rock, as well as to validate the use of the Kozeny equation in combination with a Breakage law for cohesionless materials. The latter rock is instead used with the goal to assess the role of cement disaggregation on the growth of surface area and the consequent reduction of the rock permeability. It is shown that for such materials, the mechanical response, grain size distribution (gsd) evolution and permeability reduction can be predicted satisfactorily for different loading paths by linking Breakage Mechanics formulation to a standard permeability equation for granular media, i.e., the Kozeny equation, which here has been expressed as a function of the Breakage internal variables through harmonic mean grain size. In particular, it has been shown that the most accurate predictions of loss of hydraulic conductivity are obtained by the permeability model relying on both porosity changes and surface area growth due to grain crushing. By contrast, it has been found that permeability models relying only on porosity changes tend to severely overestimate the value of rock permeability (thus providing an upper bound), while those accounting for the effect of surface area change due to cement disaggregation underestimate its values (thus representing a lower bound for the expected permeability of a comminuted rock). The performance of the proposed hydro-mechanical models suggests that this approach may provide enhanced capabilities for the simulation of evolving porosity and permeability in reservoir rocks, thus representing a convenient tool to connect the large-scale implications of fluid depletion to their inherent microscopic causes. As a result, this work motivates future analyses that can further expand the scope of microscopic processes taken into account for the computation of permeability reduction, such as the rupture of asperities and the generation of fines prior to macroscopic yielding, the heterogeneity of the mineral content in the cement and skeleton phases and the possible emergence of heterogeneous deformation processes in the form of compaction and/or shear bands.
p r =1 kPa is the reference pressure, K is the non-dimensional material parameter representing the bulk stiffness for grains or cement bond, G is the non-dimensional material parameter representing the shear stiffness for grains or cement bond and A = K ε e v /2 + 1. Depending on the material's elastic response, either linear, nonlinear or even hybrid forms (i.e. nonlinear function for the volumetric part and linear function for the shear part and vice versa) of the reference energy can be selected. In this paper, nonlinear elasticity had been used for cohesionless materials for both volumetric and shear reference energy functions while a combination of nonlinear volumetric and linear shear functions has been used for the cemented granular rocks.
In addition, the flow rules for this model derived from the dissipation potential [22] are listed in Equation (20).Ḃ = 2λ ( whereλ is a non-negative plastic multiplier, ε p v and ε p s are the plastic volumetric and shear strains, respectively, and ω B is a model parameter coupling grain crushing and the volumetric plastic dissipation associated with the reorganization of crushed fragments [35].