Shear-enhanced compaction analysis of the Vaca Muerta formation

Laboratory measurements on Vaca Muerta formation samples show stress-dependent elastic behavior and compaction at representative in-situ conditions. Experimental results show that the analyzed samples exhibit elasto-plastic deformation and shear-enhanced compaction as the main plasticity mechanism. These experimental observations conflict with the anticipated linear-elastic response prior to the brittle failure reported in several works on the geomechanical characterization of the Vaca Muerta formation. Therefore, we present a complete laboratory analysis of samples from the Vaca Muerta formation showing experimental evidence of nonlinear elastic and unrecoverable shear-enhanced compaction. We also calibrate an elastoplastic constitutive model using these experimental observations; the resulting model reproduces the observed phenomena adequately.


Introduction
Reservoir rocks' nonlinear and heterogeneous nature is typically simplified when analyzing deformation and failure during oil-and-gas well operations such as drilling, hydraulic fracturing, and production [8].Consequently, geomechanical engineers routinely assume that the material response is linear-elastic until reaching brittle failure; this assumption allows them to use straightforward analytical approximations to model stress distribution around a wellbore [10].These oversimplifications are widely used to model wellbore stability problems and hydraulic fracture [16].Typically, the linear elasticity assumption in unconventional reservoirs has been justified based on their brittleness [15], implying that reservoir rocks favorable to hydraulic fracturing treatments could be accurately modeled using linear elastic fracture mechanics [11].This practice generally produces acceptable results when the reservoir rock exhibits a strong linear elastic response.
The Vaca Muerta formation is Argentina's main unconventional reservoir in the Neuquén Basin.This formation is composed of sedimentary rock rich in organic mudstones, limestones, and marls; it was deposited in a distal ramp [23].The mineralogy of this type of reservoir rock is dominated by calcite, quartz, mica, pyrite, and clays.Typically, Vaca Muerta mudrock is characterized as a linear elastic material [27] in rate-independent mechanical problems or as a visco-elastic material [12] in geological time-dependent basin modeling problems.However, during routine triaxial testing to characterize the elastic properties of a set of Vaca Muerta samples, we observe shear-induced compaction; this compaction as a plasticity-driven mechanism may explain specific field observations such as inefficient fracture initiation, unexpected wellbore production underperformance due to unexpected hydraulic fracture geometry, or proppant placement [20,29].
The compaction of porous rocks is often explained as the closure of porosity due to increasing effective stress, assuming that solid constituents have negligible compressibility.This mechanism is a factor in the life-cycle of conventional reservoirs in the form of permeability reduction [3] and subsidence [9].In addition, the regional in-situ stress state impacts the rock's failure mode as demonstrated in numerous studies found in the literature addressing the brittle-ductile transition [2,19,[30][31][32].During laboratory testing, as confining pressure increases, the rock failure evolves from void volume creation through the formation and opening of micro-cracks that finally coalesce into a macroscopic material discontinuity known as brittle faulting [13].This failure mode involves the accumulation of irreversible plastic volumetric strain from grains rearrangement, the collapse of the pore volume or grain fragmentation and it is known as shear enhanced compaction [6].Plastic deformation in mudrocks is typically attributed to the abundance of organic matter and clay inside the reservoir rock matrix [28].However, the laboratory experiments conducted in this work in the Vaca Muerta mudrock samples suggest that shear-enhanced compaction, possibly controlled by grain displacement, is also a predominant mechanism of unrecoverable volumetric plastic deformation.Here, we characterize the compaction mechanism experimentally, conducting a series of triaxial tests over a set of samples from the Vaca Muerta formation.
We adopt a phenomenological nonlinear elastoplastic constitutive model that adequately captures the observed material response; we calibrate its constitutive parameters for use in engineering applications.From all available phenomenological constitutive models, we select the Modified Cam-Clay (MCC) model [21,22] as it captures the main experimental observations: nonlinear elasticity and isotropic hardening evolution of the yield surface, which translates into unrecoverable volumetric deformation.Although, MCC was originally developed for characterizing the soils' critical state, it was extended to other types of cohesive-frictional materials [4,7].In addition, we implement an implicit integration algorithm [1,17] and simulate a triaxial test response at a Gauss point.We structure the manuscript as follows: a brief introduction in this section; Section 2 gives notation and preliminary definitions that are used throughout this work; Section 3 briefly introduces the plasticity theory in the context of thermodynamics of deformable solids; Section 4 addresses the continuous elastoplastic formulation, Section 5 discusses experimental evidence of elastoplastic response on Vaca Muerta mudstone, Section 6 details the numerical implementation, and, finally, Section 7 discusses our conclusions.

Notation and preliminary definitions
We denote second-order tensors (i.e., matrix arrays that satisfy certain change of basis rules) using bold Greek symbols or letters, fourth-order tensors with upper case blackboard-bold letters (e.g., C)), vectors (i.e., first-order tensors) with lower case italic, bold letters, and scalars with lower case Greek symbols or letters.In addition, we adopt Einstein's summation convention (given i = 1, 2, 3; a i b i = a 1 b 1 + a 2 b 2 + a 3 b 3 ) to perform componentwise tensor operations such as contractions (vector inner products) or double contractions (second-order-tensor inner products).The symbol ⊗ denotes the tensor product (outer product).We assume that all the operations are performed in Euclidean space (zero curvature space).Thus, we make no distinction between covariant and contravariant basis vectors.Unit basis vectors in the Euclidean Space R 3 are denoted as e i , i = 1, 2, 3 with the property e i • e j = δ i j and e i • e i = δ i i = 3 for i , j = 1, 2, 3, where δ i j ∈ R 3×3 is the Kronecker Delta, where δ i i = 1 and δ i j = 0 if i ̸ = j .Therefore, we write first-order, second-order, and fourth-order tensors as follows,

Remark 1 (Index expansion).
In what follows, we work in a three-dimensional Euclidean Space.Thus, the range for indexes will be omitted in the rest of this manuscript.
The transpose of a second-order tensor σ ∈ R 3×3 using index notation is defined as, and the trace of a second-order tensor σ is given by tr(σ) = σ : 1, where 1 = δ i j e i ⊗ e j .Given two vectors a, b ∈ R 3 and two second-order tensors σ, ε ∈ R 3×3 , the contraction and the double contraction are, We often make use of the symmetric fourth-order unit tensor defined as, The symmetric fourth-order unit tensor admits the following decomposition: where I v and I d are the volumetric and deviatoric contributions given by, The operator ∥ * ∥ : * → R denotes the Frobenius norm for either vectors or second-order tensors.For vectors, ∥a∥ = a • a = a i a i , whereas for second-order tensors ∥σ∥ = σ : σ = σ i j σ i j .We usually work with variables that evolve during a pseudo-time increment ∆t = [t n+1 , t n ], ∀n ∈ N + .The derivative of a tensor or scalar field with respect to the pseudo-time t is denoted as: We adopt the rock mechanics convention in which compressive stresses are positive and tensile stresses negative.
In addition, we denote the Effective Cauchy's Stress Tensor by σ and the Infinitesimal Strain Tensor by ε.The Hydrostatic Stress is p = 1 3 tr(σ) = σ : 1 and the Deviatoric Effective Stress s = σ − p 1 with the property that tr(s) = 0. Additionally, for a symmetric second-order tensor σ, we define the following tensor invariants, The stress invariants for the deviatoric symmetric second-order tensor are, Finally, the deviatoric stress invariants relate to the stress invariants as follows,

Thermodynamics of deformable continua
This section introduces the mathematical theory of plasticity in the context of the thermomechanics of deformable bodies (see [18]).Let B be a continuum body with volume V and boundary ∂B.Let us define on B the following scalar fields: ρ, θ, u, s, and R representing the mass density, temperature, internal energy, entropy, and heat density, respectively.We assume that small deformations properly describe the system's kinematics.Thus, the first (energy conservation) and second (entropy production imbalance in the form of the Clasius-Duhem inequality) thermodynamic laws in their local form adopt the following differential expressions, 1.Energy conservation: where σ, ε ∈ R 3×3 are Cauchy's stress and strain rate tensors, respectively, ∇ = ∂ ∂x i e i is the gradient operator, and q is the heat flux vector acting on ∂B.

Entropy production imbalance:
Substituting ( 12) into (13) we get the Clausius-Duhem Inequality, We define the Helmholtz free energy per unit of mass as, Recalling that (by distributive property of the ∇ operator over vector and scalar fields), and the change of variables g = ∇θ, we rewrite (14) as, Assumption 1 (Thermodynamic evolution).We assume that the system's evolution is isothermal, leading to a purely mechanical formulation for (16), Remark 2. Every mechanical system must satisfy the Clausius-Duhem inequality of (17) where the system's energy in thermodynamic equilibrium is characterized by a series set of state variables.
Axiom 1 (Local-state postulate).Given an arbitrary system evolution, the energetic state of a continuum medium is characterized by the same state variables that are fixed at the equilibrium state, and it is rate independent.Therefore, the Helmholtz free energy (15) is given by, where η = η k , ∀k ∈ N + is an internal variable set describing the material's dissipation and its free energy rate is, The symbol " * " denotes the contraction compatible with ∂Ψ ∂η and η.

Reversible and irreversible thermodynamic processes: elastoplasticity
We characterize the deformation evolution of a material as reversible (non-dissipative) and irreversible (dissipative) by formalizing these concepts using the following definitions: Definition 1 (Reversibility: elasticity).A deformation process is non-dissipative, reversible or elastic if and only if the internal variables remain constant throughout the deformation process (i.e., η = 0).Thus, (17) Letting Ψ = W , where W is the strain energy density, we express the hyperelastic constitutive models as We define dissipative thermodynamic processes by adopting an additive decomposition of the strain tensor, Assumption 2 (Strain-tensor additive decomposition).Let ε ∈ R 3×3 be the strain tensor at a material point x ∈ B, which admits the following decomposition, where ε e and ε p are the elastic and plastic components of the strain tensor, respectively.

Definition 2 (Dissipative Processes).
A thermodynamic deformation process is dissipative, irreversible, or plastic if the free energy u can be decomposed in a strain energy density W and a latent energy density V .Thus, defining Therefore, the Clausius-Duhem inequality (17) becomes, Remark 3 (Reversible and irreversible processes).A system's evolution is reversible if and only if there exists an isomorphism between the initial and the final state; otherwise, the evolution of a system is irreversible.
For irreversible processes, the Clausius-Duhem inequality is satisfied for more than one set of state variables.Thus, we enforce uniqueness using the Maximum Plastic Dissipation Condition.
Let E σ be a closed and convex set defining the admissible state variables E σ given by where F f : R 3×3 ×R n → R bounds the admissible stress states (yield function).This assumption defines a unique state variable set σ, η, εp such that, Remark 4 (Elastic Domain and Flow Surface).The set of admissible states E σ admits a partition and ∂E σ is the flow surface defined as Proposition 1 finds the state variables that maximize D p subject to the constraint F f (σ) = 0 (Consistency Condition) and the complementary Kuhn-Tucker conditions [26].We solve this optimization problem by introducing a Lagrangian (cost function) that transforms the maximization into a minimization problem (i.e., let −D p (τ, κ, εp )) where L : R 3×3 × R n × R + → R is the Lagrangian and γ ∈ R + is the Lagrange multiplier.Henceforth, satisfying the Maximum Plastic Dissipation Condition is equivalent to solving the following minimization problem: The necessary optimal condition is ∇L | (σ, η) = 0 and the Lagrangian is convex (concave); thus, we can deduce, ∂L ∂σ ∂L ∂η where H is the Hardening Modulus given by, From ( 28) and ( 29) and considering the following change of variable γ = λ, the Generalized Associative Flow Rule and the Generalized Hardening Law are [24], • Generalized Associative Flow Rule: • Generalized Hardening Law: Remark 5 (Non-Associative Flow Rules).Proposition 1 is a sufficient condition that is too restrictive in general, but it inherently induces a generalized hardening law and a generalized associative flow rule.In practice, the associative plastic flow rule appropriately characterizes materials with crystalline micro-structural composition.Although associative flow rules are used to characterize the plastic flow of granular materials, non-associative flow rules could be more suitable.Therefore, the Clausius-Duhem Inequality in its local form should be verified independently when using non-associated plastic flow rules.Typically, a non-associated flow rule is defined as These heuristic definitions seek to reconcile experimental observations with simulations by relaxing the overly restrictive Proposition 1, while their major weakness is their ad-hoc nature.

Continuous elastoplastic constitutive model
Elastoplastic constitutive models for materials should capture the following experimental observations [24]: 1.For loads below a threshold that defines a flow criterion, the material response is reversible (elastic).
2. Once the material reaches the limit condition, the deformation becomes partly irrecoverable (plastic).
3. Plastic deformations evolve the failure state as described by a hardening law.4.During unloading, the material response is elastic.5.The material response during the whole deformation process is quasi-stationary.
6.The material is thermodynamically stable.
We describe the material's elastic response by adopting a linear Hookean model given by, where C e is the fourth-order elasticity tensor defined by the bulk modulus K , and the shear modulus G as, with K and G adopt the following expressions [17], where κ is the volumetric deformation recovery, φ is the rock's total porosity, and ν is Poisson's ratio.
Remark 6 (Pressure dependency of K and G).During laboratory testing conducted in Vaca Muerta samples, we observed pressure dependence on the bulk K and shear G moduli.We capture that phenomena including isotropic part of the effective stress tensor p in (36) and (37).However, some porous rocks in specific stress state ranges could be modelled as linear and non-pressure dependent.
Remark 7 (Coupling K and G).The coupling of elastic shear and volumetric moduli may lead to energy dissipation under cyclic loading.However, non-conservation is not an issue for monotonic loading.In addition, the definitions for K and G in (36) and (37) are widely used in practice; thus, we adopt these expressions.
We also assume that the porosity evolution during loading satisfies the following state equation in rate form, where ψ represents the porosity degradation rate during volumetric deformation (see 12).The mathematical description of an elastoplastic constitutive model involves three main components [14]: 1. Yield Function: describes the location of points where the material develops irrecoverable deformation.
2. Flow Rule: characterizes the evolution of irrecoverable deformations 3. Hardening Law: typifies the evolution of the yield function throughout the plastic evolution.

Yield Function: Modified Cam-Clay Model
Cam-Clay models [21,22] are widely used for plasticity characterizations of the stress-strain response of cohesivefrictional materials subject to three-dimensional stress states [4,7].These simple models can realistically represent the compaction and dilation responses of porous materials [1,17].Cam-Clay models capture the typical pressure sensitivity and hardening of cohesive-frictional materials, requiring few parameters that standard laboratory testing procedures can characterize (see [21,22]).Recently, experimental data on limestone rocks that exhibit compaction and dilation during laboratory triaxial testing was reported in [25].The same material response was observed in Vaca Muerta mudstone rock during triaxial testing reported in this work.Thus, we select a constitutive model that captures this nonlinear response for practical geomechanical applications; see Figure 1 for a schematic representation of the MCC yield criterion.
We adopt a Modified Cam-Clay Model defined in the p − q plane as, where p and q represent the hydrostatic pressure and the deviatoric part of the effective stress tensor given as, and M , p c define the Critical State Line (CSL) slope and the hardening parameter, respectively.Alternatively, we express (39) in terms of the stress invariants as, where Figures 1(a) and (b) represent (39) and (42), respectively.
Remark 8 (Hardening parameter).The state variable p c is a hardening parameter that describes the yield surface's evolution as the hydrostatic pressure grows.

Flow Rule
We adopt an associative flow rule, which expresses the rate of plastic deformation as, where the flow direction is given by (see Appendix A.1), in which n = s ∥s∥ and λ > 0 is the plastic multiplier.
Remark 9 (Non-associativity).Although associative flow rules are used to characterize the plastic flow of granular materials, non-associative flow rules on the volumetric component of εp could be more suitable for certain type of rocks.Nevertheless, assuming associative flow to capture the plastic flow in Vaca Muerta reproduce the experimental observation without adding any additional heuristics.

Hardening Law
We define the hardening law in rate form as, where εp v = εp : 1 is the volumetric component of the plastic strain rate tensor and χ is given by In ( 46   Triaxial compression tests for vertically oriented samples were conducted at W. D. Von Gonten Engineering laboratory facilities.Each of the samples was tested in a modern servo-controlled triaxial loading system in accordance with ASTM testing standards.Laboratory tests were performed under room temperature and drained conditions using Teflon isolating jacket, cantilever radial displacement gauge, and four equally spaced axial displacement transducers (LVDT).Before conducting the triaxial test, three hydrostatic cycles and uniaxial strain compression cycles were performed before the triaxial test (see Figure 2).The initial cycles of hydrostatic compression consolidate each sample and minimize the effect of coring-induced microcracks and other artifacts.The uniaxial-strain test seeks to evaluate the evolution of elastic properties as a function of the confining pressure.Before the failure, a triaxial compression test was used to investigate the elastic properties and the post-elastic stress-strain behavior, including dilatant and compactive plasticity.We test four Vaca Muerta mudstone samples, and Table 1 summarizes their properties.

Material parameters estimation from laboratory tests
We assume that the total volume of the porous medium admits an additive decomposition V = V v +V s (see [5]) where V v and V s are the pore and solid volumes, respectively (see Figure 3); therefore, V v and V s are,

Assumption 3 (Matrix volume conservation). At any load step, the matrix volume remains constant; thus,
This assumption is relevant as we are not destroying solid volume at laboratory conditions.

Material response during hydrostatic cycling
During the hydrostatic cycling stage, deviatoric stress q = σ 1 − σ 3 remains constant and hydrostatic pressure p = 1 3 (σ 1 + 2 σ 3 ) is monotonically increased from 200 psi to 2600 psi.From this test, we determine the initial volumetric modulus during unloading K considering the slope of the unloading curve from p against the ε v chart.
Additionally, κ and γ material parameters are determined by applying Assumption 3 to calculate φ and analyzing its variation as a function of p and identifying the slopes of the loading and unloading curves.The initial hardening parameter p 0 c is determined from the intersection of the loading and unloading curves in the φ − p plot.Remark 11 (Hysteresis during unloading).During unloading/reloading cycles, an insignificant amount of hystere-sis develops.We disregard this behavior for parameter interpretation and consider only the first unloading response to define each material parameter and the initial hardening p 0 c .
Remark 12 (Maximum hydrostatic pressure).The maximum hydrostatic pressure during cycling should be the confinement pressure at triaxial conditions.We subject each sample to 2500 psi of peak hydrostatic pressure for the data we present.Otherwise, the triaxial stress path may induce unwanted isotropic compression with a certain amount of plastic compaction.Therefore, the sample might show a different response at failure.

Material response during triaxial test
A triaxial test was conducted after the hydrostatic and uniaxial stages until failure.Each sample was subjected to a confinement pressure (σ c ) of 2500 psi during the test.Figures 8 to 11 show Vaca Muerta mudstone material response during the triaxial test.These figures show a linear elastic response until 0.05 % of axial deformation.Within the linear elastic range, we estimate Young's Modulus and Poisson's Ratio; we summarize these elastic constants in Table 3.
Parameter 2-44v 3-14-2v    Finally, the material dilates until localized shear failure.Nonlinear compaction is evident when we plot the hydrostatic stress p against ε v where the graph departs from the theoretical hydrostatic linear trend (point p c in p vs ε v charts).Thus, plastic volumetric strains develop due to compaction until a limit point where the material dilates at constant hydrostatic pressure.We observe a similar response when representing the φ evolution, when considering Assumption 3 with increasing hydrostatic pressure.Initially, the volumetric plastic deformation induces a monotonic porosity reduction until a critical point where the material dilates and, ultimately, fails in shear.The state equation ( 38) describes the porosity evolution in rate form.The porosity degradation parameter ψ is conventionally derived from plotting φ against ε v .Figure 12 shows a typical porosity degradation profile of our  Clay evolution through the experimental stress path for the Vaca Muerta mudstone samples.The slope of the stress path is 3, the theoretical slope for triaxial conditions.The samples are initially loaded to an initial state p 0 , q 0 inside the initial yield locus defined by the major axis p 0 c (the initial hardening parameter captured by hydrostatic  The slope of the CSL (M ) is defined by the intersection between the triaxial test stress path and the hydrostatic pressure p at which the sample starts to increase volumetric strain at constant pressure looking at the ε v vs p chart.Thus, the CSL is a straight line from the origin to the latter intersection on the q −p projection.Table 4 summarizes the MCC constitutive parameters for the Vaca Muerta mudstone samples we tested.

Closest Point Projection Mapping for the Modified Cam-Clay Model
Following [1,14,17], let Ω ∈ R d , d = 2, 3. Consider the discretization Ω h ⊂ Ω.Let us take an arbitrary Gauss point on a finite element e ∈ Ω h .Let t ∈ [0, T ], T > 0 be a pseudo-time anf n, k ∈ N + be the pseudo-time increment and iteration counters, respectively.The incremental strain tensor ∆ε k n+1 is, Additionally, consider the trial state defined by freezing all the internal variables as: where σ n and ε n are the converged effective stress and strain tensors of the previous pseudo-time step n.The return mapping tensor equations in their general form are: Integrating ( 31) within [t n , t n+1 ] leads to the following discrete plastic-strain increment ∆ε p , where ∆λ is the discrete consistency parameter.Consider the volumetric part of σ k n+1 : where ∆ε Combining ( 56) and (57), results in the following system of scalar equations on ∆λ: Exact integration of the hardening law (45), gives: The consistency parameter ∆λ in (58))-( 59) is calculated by imposing the consistency condition on F f (∆λ): Since (59) couples the variables p k n+1 and p c , F f (∆λ) cannot be evaluated explicitly.Therefore, p k n+1 and p c are updated iteratively by a local Newton-Raphson iteration.Thus, we rewrite the first equation in (58) as, In Section 4, we define K and G as dependent on p and φ.Therefore, K and G are state variables updated at each strain increment during loading.We use an explicit integration scheme to update φ, K and G as follows: Thus, combining (49) to (65) and the derivatives of F f (∆λ) and H f (∆λ) with respect to the consistency parameter ∆λ (see Appendix A.3) to fully define Newton-Raphson iteration scheme, we obtain the following closest point projection algorithm (see [14,17]) for updating σ k n+1 , ε p n+1 and p c (see Figure 14 for a sketch of the closest point projection algorithm detailed in Algorithm 1):

Numerical Results: Triaxial test simulation
Next, we evaluate the performance of the closest point projection algorithm to reproduce a triaxial test by focusing on a single Gauss point P G and inducing a triaxial stress state.We consider a confinement pressure σ c applied to the sample's sides; at the sample's top, we apply a strain increment ∆ε, which induces a deviatoric strain increment at the Gauss point.Additionally, we fix the sample's bottom during the numerical experiment.Figure 15 shows the general setting for the numerical experiment.
We consider an initial hydrostatic stress state at P G of the form, Additionally, at P G we apply the following deviatoric strain tensor, 3. Evaluate the Cam-Clay Yield function at the trial state:
We consider material parameters within the range we obtain for the Vaca Muerta mudstone samples.Table 5 summarizes the input data we use to simulate the triaxial test stress response.Figure 16 shows the triaxial test simulation using the MCC Yield function and Closest Point Projection algorithm to update state variables at the Gauss point P G and the triaxial response for Vaca Muerta mudstone samples.The calibrated constitutive model reproduces the main features we observe in the stress-strain material response in the lab tests.This model prop-  erly captures plasticity due to shear-enhanced compaction despite using an associative flow rule for updating state variables.However, the model parametrization cannot represent the dilation response, as our state equation we introduce does not capture this phenomenon since the hardening law induces compaction as the yield function evolves.Moreover, this model predicts perfect plasticity instead of dilation after sample reaches its critical deformation.Even though some samples showed dilation and softening, these phenomena happens at the end of the triaxial test just before the shear failure of the sample.Thus, our model properly captures the main plastic response due to compaction.

Conclusions and future work
We calibrate a constitutive model capable of capturing the volumetric plasticity of the Vaca Muerta mudstone using a standard laboratory testing program.The integration algorithm inspired by [17] can properly update the state variables σ k n+1 and ε p n+1 ; this algorithm uses a closest point projection strategy proposed primarily by [14] with an associative flow rule and a simple hardening law.Although we adopt an associative flow rule, the model is capable of capturing compaction without over-predicting volumetric strains.Our mathematical model does not seek to describe the dilatant behavior of this rock as the plastic load evolves and the algorithm updates the state variables.Therefore, the material behaves as perfectly plastic after a critical compaction value.The model Proof.Since p = 1 3 σ : 1 and using index notation, the derivative with respect to the effective stress tensor is, The derivative of q with respect to the effective stress tensor is collected in the following results: Lemma 2. Let σ, s ∈ R 3×3 be a second-order tensor and its deviatoric part.Let σ i j and s i j , i , j = 1, 2, 3 the components of σ and s respectively.The following holds, By Proof.We expand the left-hand side of the equality using index notation, the definition of the norm, and double contraction for second-order tensors,  Thus, the derivative of F f with respect to the effective stress tensor adopts the following expression, Proof.The result follows by applying the chain rule for derivation,

Figure 1 :
Figure 1: Cam-Clay Yield Surface represented in different spaces.
) γ, κ are material parameters typically obtained in hydrostatic cycling tests.Remark 10 (On the γ and κ definitions).γ and κ are material parameters representing compaction and bulk volumetric recovery measured in the hydrostatic cycling test.The required number of hydrostatic cycles applied during testing induces a linear trend from loading and unloading cycles.At this point, the sample is compacted at a certain level where artifacts from sample extraction for the core are mitigated.

Figure 2 :
Figure 2: General laboratory test program for mudstones and limestones.

Figure 3 :
Figure 3: Additive decomposition of a porous medium volume.

Figure 12 :
Figure 12: Typical Porosity degradation profile at increasing volumetric strain.

Figure 14 :
Figure 14: Internal variables update by the Closest Point Projection Algorithm

Figure 15 :
Figure 15: Triaxial test general setting for numerical experiments.

= δ i k δ l j − 1 3 δ i j δ kl . Proposition 3 .2 ∥s∥ = 3 2 s
Lemma 1 and Proposition 2 and expanding the left-hand side of the previous equality, Let q = 3 : s be a measure for the deviatoric stress tensor; thus, the following holds,

= 3 2 s pq s pq − 1 2 σ i j − 1 3 σ l l δ i j = 3 2 s pq s pq − 1 2pq s pq − 1 2Proposition 4 (
σ i j − p δ i j ) and contracting indexes, we get the desired expression.σ i j − p δ i j e i ⊗ e j collect the derivatives of the modified Cam-Clay yield function is collected in the following, Cam-Clay Yield Function Derivative).Let F f (σ) : R 3×3 → R be the Cam-Clay yield function then, F f (σ) := q 2 M 2 + p (p − p c ), ∀M , p c > 0.

Table 1 :
Vaca Muerta sample properties as received.

Table 2 :
Material parameters for each Vaca Muerta sample.
Figures 4 to 7 show the material response during hydrostatic cycling for Vaca Muerta samples and the resulting parameters K , κ, γ and p 0 c .Table2summarizes the interpreted material parameters from the hydrostatic test.

Table 3 :
Elastic Constants for each Vaca Muerta sample.

Table 5 :
Numerical Experiment.Constitutive and Material Parameters