High‐Resolution Mapping of Yield Curve Shape and Evolution for High‐Porosity Sandstone

Understanding the onset and nature of inelastic deformation in porous rock is important for a range of geological and geotechnical problems. In particular for sandstones and siliciclastic sediments, which often act as hydrocarbon reservoirs, inelastic strain can significantly alter the permeability affecting productivity or storativity. The onset of inelastic strain is defined by a yield curve plotted in effective mean stress (P) versus differential stress (Q) space. Yield curves for porous sandstone typically have a broadly elliptical shape, with the low‐pressure side associated with localized brittle faulting (dilation) and the high‐pressure side with distributed ductile deformation (compaction). However, recent works have shown that, for different porous rocks, the curve shape can evolve significantly with the accumulation of inelastic strain. Here yield curve shape and evolution of two high‐porosity sandstones (36–38%) is mapped along different loading paths using a high‐resolution technique on single samples. The data reveal yield curves with a relatively shallow geometry and with a compactive side that is partly comprised of a near‐vertical limb. Yield curve evolution is found to be strongly dependent on the nature of inelastic strain with samples compacted under a deviatoric load (i.e., with a component of shear strain) having peak stress values that are approximately 3 times greater than similar porosity samples compacted under a hydrostatic load (i.e., purely volumetric strain). These results have important implications for predicting how the strength of porous rock evolves along different stress paths, which differ in reservoirs during burial, fluid extraction, or injection.


Introduction
There is widespread interest in the cause of porosity loss from sandstones and siliciclastic sediments, primarily because of the impact on fluid flow and permeability reduction in reservoirs. During burial, sandstones predominantly lose their primary porosity by mechanical and/or chemical compaction. Although chemical compaction, particularly through pressure solution and cementation, may significantly reduce fluid flow in deep reservoirs, the initial stages of burial are typically dominated by mechanical processes such as grain rearrangement, grain crushing, and localized faulting (e.g., Worden & Burley, 2003). Understanding the nature of mechanical deformation in porous rock is important for both geological and geotechnical problems including faulting and deformation band formation (Antonellini et al., 1994;Aydin & Johnson, 1978;Ballas et al., 2013;Griffiths et al., 2016), borehole stability (Cuss et al., 2003a;Dresen et al., 2010), and reservoir subsidence and compaction (Fisher et al., 1999;Makowitz & Milliken, 2003;Nagel, 2001;Pijnenburg et al., 2019).
The onset of mechanical deformation in porous rock is typically defined using a yield curve plotted in P-Q space (Figure 1a), where P is the effective mean stress P ¼ σ1þσ2þσ3 3 −P f À Á and Q is the differential stress (Q = σ 1 − σ 3 ). Here P f is the pore fluid pressure and σ 1 , σ 2 , and σ 3 are the principal stresses. The yield curve separates regions of elastic deformation (i.e., recoverable) from regions of inelastic deformation (i.e., permanent). The response of a porous rock at the onset of yield can be broadly separated into two regimes. At low effective pressures, inelastic deformation is associated with dilatancy and localized faulting (e.g., Menéndez et al., 1996). Under higher effective pressures the deformation transitions into a regime of cataclastic flow and distributed compaction (e.g., Curran & Carroll, 1979). This transition from localized to distributed deformation is often referred to as the low temperature brittle to ductile transition (e.g., Rutter & Hadizadeh, 1991;. Many of the principles of porous rock deformation are based on those of critical state soil mechanics (e.g., Schofield & Wroth, 1968;Wood, 1991). However, rocks are not soils, and thus, the micromechanisms of deformation will be different, which must be considered when applying these principles. Traditionally, yield curves of porous rock have been mapped by performing a suite of axisymmetric compression tests (σ 1 > σ 2 = σ 3 ), where multiple samples of a given rock type are axially loaded under different effective pressures. This has been done on a variety of different rock types including sandstone (Baud et al., , 2006Baud, Zhu, et al., 2000;Cuss et al., 2003b;Louis et al., 2009;Wong et al., 1997), limestone Baud, Schubnel, et al., 2000;Cilona et al., 2014;Vajdova et al., 2004), volcanic rocks (Zhu et al., 2011), and dehydrated serpentinite (Rutter et al., 2009). Most yield curves have a broadly elliptical shape with the lowpressure side of the curve associated with dilatancy and the high-pressure side with compaction ( Figure 1a). As yield of porous rock can be achieved under purely hydrostatic stress, the yield curve intersects the P axis at the pressure where hydrostatic yield occurs. This point (P*) marks the onset of pore collapse and grain crushing (e.g., . Traditional investigations, using axisymmetric compression tests on multiple rock samples, have provided many insights in to the general shape of yield curves and the nature of deformation in porous rock. However, some studies have found that the yielding behavior of sandstones Fortin et al., 2005) and porous volcanic rocks (Heap et al., 2015) can be better approximated using linear yield curves rather than ellipses. Alternative models to the traditional elliptical yield cap model of critical state soil mechanics have been used to describe these linear yield curves (Guéguen & Fortin, 2013;Zhu et al., 2010). More recently, new methodologies have been developed to try and map precisely the shape of yield curves for single samples of porous rock in order to reduce the effects of sample variability and better constrain their geometry (Bedford et al., 2018;Tembe et al., 2007). These studies found yield curves to be elliptical in shape but with a near-vertical region close to the proximity of P*, which was termed by Bedford et al. (2018) to be a region of shear-insensitive compaction (Figure 1a). The variations in yield curve geometry observed in previous investigations (e.g., elliptical, linear, and vertical limbs) highlight the importance of gathering more high-resolution data to constrain better the yield curve shape for a variety of different porous rock types.
As well as a yield curve that identifies the stress conditions at which the onset of inelastic strain will occur, a comprehensive description of porous rock deformation also requires a work hardening rule, which describes how the yield function changes as inelastic strain is accumulated (Paterson & Wong, 2005). During Figure 1. A schematic illustration of yield curves plotted in P-Q space. (a) Yield curves for porous sandstone are broadly elliptical in shape with the low-pressure side associated with dilatancy and the high-pressure side with compaction. The point where the curve intersects the effective mean stress axis is referred to as P*. Recent works have shown that region of the yield curve in the vicinity of P* is near vertical (Bedford et al., 2018;Tembe et al., 2007) and associated with shear-insensitive compaction. (b) A family of yield curves can be defined for different porosities. The peak of these curves is joined by the critical state line to separate the regions of dilation from compaction. (c) A 3-D representation of the yield envelope where the third axis is porosity (Φ) or sometimes porosity multiplied by grain size. The yield curves space out along the porosity axis, and their respective P* points are joined by the normal consolidation line. deformation, the porosity will either increase or decrease as a result of compaction or dilation causing the yield curve to expand or contract. A family of yield curves can therefore be defined to represent all the subsequent changes in porosity during deformation ( Figure 1b). Traditionally, the crests of these curves are joined by the critical state line (CSL) that separates regions of dilation from compaction. The family of curves can also be visualized as a three-dimensional surface, where the third axis is porosity (or porosity multiplied by the grain radius) and the surface is formed by the continuum of individual yield curves that space out along this axis (Figure 1c). The normal consolidation line, also referred to as the hydrostat, represents the compaction behavior of a porous material to purely hydrostatic loading. In soil mechanics, isotropic hardening is assumed, which produces a successive family of yield curves that symmetrically expand with compaction (Schofield & Wroth, 1968). Therefore, the typical view is that this family of curves is unique with a universal shape that does not evolve with deformation. However, experimental studies have begun to show that this concept does not apply to porous rock deformation and that significant yield curve shape evolution can occur with the accumulation of inelastic strain (Baud et al., 2006;Bedford et al., 2018). Plasticity models, which build on the principles of soil mechanics but allow an evolving yield criterion, have been developed to account for the effects of inelastic volumetric strain (e.g., Carroll, 1991;DiMaggio & Sandler, 1971) and inelastic shear strain (Grueschow & Rudnicki, 2005). However, more empirical data are needed to test the effects of different types of strain on yield curve evolution and to constrain the effects of stress path on porous rock strength.
The study of Bedford et al. (2018) is, to the best of our knowledge, the first to separate the effects of inelastic volumetric strain and inelastic shear strain on yield curve evolution. This was done on porous bassanite, a material produced from the dehydration of gypsum and comprised of relatively weak platy grains. A stress-probing methodology on single samples was used to perform a series of overconsolidation tests where samples are deformed beyond their elastic limit under different hydrostatic and deviatoric loading paths. Results showed that a sample compacted under a deviatoric load (i.e., with a component of volumetric and shear strain) had a yield curve with a peak that is almost double that of a sample that had been hydrostatically compacted (i.e., purely volumetric strain) to the same porosity. The difference in yield curve shape is qualitatively attributed to the heterogeneous microstructure and formation of multiple sets of conjugate shear bands that developed during deviatoric loading, highlighting the importance of understanding microstructural controls on rock strength. Bedford et al. (2018) also document a shape change with the formation of a plateau at the peak of the yield curve during deviatoric loading. They describe this region as a critical state zone rather than the CSL that is typically considered in porous rock deformation (e.g., Figure 1b). After the work by Bedford et al. (2018) on porous bassanite, it is apparent that the effects of yield curve evolution need to be studied more extensively for other porous rock types, where the constitutive behavior has been constrained previously, to see whether the same evolutionary behavior along different stress paths is also applicable. Since the seminal work of Wong et al. (1997), there have been many investigations into the deformation behavior of porous sandstone (e.g., Wong & Baud, 2012, and references therein). In this study, we therefore apply the stress probing methodology of Bedford et al. (2018) to test if the yield behavior and subsequent evolution observed for porous bassanite is also observed for sandstone. We map in high resolution the yield curves of two high-porosity sandstones (36-38%). The postyield evolution of the curves is investigated in response to both hydrostatic and deviatoric loading and the associated microstructures from these different loading histories are analyzed.

Experimental Materials
The two porous sandstones used in this study are Boise and Idaho Gray sandstone. Blocks of both rock types were purchased from Kocurek Industries Inc. (USA). A petrophysical description of the rocks is provided in Table 1. In order to investigate postyield evolution, it is important that the sandstones are relatively weak so that they will yield well within the limits of the deformation apparatus to allow for them to be inelastically overconsolidated. Boise sandstone was chosen for this investigation because of the relatively low P* values (42-105 MPa) that have previously been reported for this rock type (Baud, Zhu, et al., 2000;Cheung et al., 2012;David et al., 1994;Wong et al., 1997;Zhang, Wong, Yanagidani, et al., 1990). The differences in the P* values recorded in previous studies are attributed to variations in porosity (25-35%) and grain size (0.28-to 0.46-mm radius) between the different blocks of Boise sandstone used in the previous investigations. In this study, initial porosities were determined using He-pycnometry, and Boise sandstone was found to have a starting porosity of 37.8% (average of 14 measured samples), which is the highest that has been reported for this rock type. Idaho Gray sandstone was chosen because of its similar petrophysical properties to Boise sandstone, which is a good indication that it should also have a relatively low P*, allowing for the postyield evolution to be investigated. There are some slight differences, in that Idaho Gray sandstone has a lower measured starting porosity of 36.2% (average of 15 measured samples), a slightly coarser grain size and is slightly more quartz-rich than Boise sandstone (see Table 1).
Cores (50-mm length × 20-mm diameter) were drilled in the same orientation from their respective blocks and were then precision ground to square the ends to tolerance (±0.01 mm). Boise sandstone displays weak bedding, and cores were drilled perpendicular to bedding. No obvious bedding or anisotropy is observed in the Idaho Gray sandstone. The porosity of each sample was measured before being placed into a 2.5-mm wall thickness polyvinyl chloride jacket prior to insertion in to the deformation apparatus.

The Deformation Apparatus
The triaxial apparatus used (supporting information Figure S1) was designed and built in the University of Liverpool and is capable of experiments under conditions of up to 250-MPa confining pressure, 200-MPa pore fluid pressure, and a 300-kN differential load (Faulkner & Armitage, 2013). The confining pressure and pore pressure are servo-controlled and can be measured to resolution of better than 0.01 MPa. The servo-controlled pore pressure pump can also be used as a volumometer to monitor changes in pore volume to a resolution better than 0.1 mm 3 . The pore fluid used in this study is argon as it is chemically inert, which was chosen to minimize the effects of any time-dependent stress corrosion cracking (e.g., Atkinson, 1984;Atkinson & Meredith, 1981;Heap et al., 2009). The pore fluid pressure was held constant at 20 MPa for all tests in this study. The axial load is generated by a servo-controlled electromechanical piston and is monitored by an internal force gauge with a resolution of better than 0.03 kN.

The Experimental Procedure
In order to characterize the postyield evolution of the different sandstones, it is first important to determine the initial yield characteristics of the starting materials. This is done in the "traditional" sense by using multiple samples to determine the yield curve. First, the onset of hydrostatic yield (P*) is constrained by incrementally increasing the effective mean stress on a sample while monitoring the relative pore volume change with the volumometer. The effective pressure was increased to >200 MPa for both rock types in order to map out the hydrostat in detail, and porosity evolution was also tracked during unloading. Then the rest of the yield curve is determined by axially loading different samples, at a constant displacement rate of 1 μm/s (strain rate = 2 × 10 −5 s −1 ), over a range of different effective pressures less than P* and monitoring the stress and pore volume evolution. Yield can be identified in these tests by the deviation from quasi-linear elastic loading (e.g., inset of Figure 2).
Once the yield characteristics of the starting material have been determined using multiple samples, the initial and postyield evolution is investigated along different stress paths using the stress-probing methodology of Bedford et al. (2018). This methodology, outlined in the subsections below, allows for the yield curve of a single sample to be mapped in high resolution, which has the distinct advantage that a complete yield curve can be mapped using very little material. This would be beneficial in cases where the amount of material available for testing is restricted, for example, recovered core from a borehole. This individual sample can then be overconsolidated along a given stress path, in order to accumulate inelastic deformation, and the new yield curve can then be mapped in the same way. Yield curve evolution is investigated along two

10.1029/2018JB016719
Journal of Geophysical Research: Solid Earth different stress paths: hydrostatic overconsolidation (isotropic loading involving only inelastic volumetric strain) and deviatoric overconsolidation (along an axisymmetric compression loading path that has a component of inelastic shear and inelastic volumetric strain). The details of each loading trajectory are outlined in the subsections below. Figure 2a shows a schematic summary of the loading path for the hydrostatic overconsolidation tests. Initially, a sample is loaded hydrostatically beyond P* so that it accumulates some inelastic volumetric strain and sits on a new yield curve (Figure 2a, step 1). As the rock has accumulated inelastic strain, it essentially becomes a new rock with a corresponding new P* value (labeled as P* I in Figure 2a). The confining pressure is then incrementally decreased so that stress state returns to the elastic regime inside this new yield curve. Between each increment of confining pressure reduction, the yield curve is probed by axially loading the sample (loading rate = 1 μm/s) until yield is achieved ( Figure 2a, step 2). As soon as yield is reached, identified by the deviation from quasi-linear elastic loading (inset of Figure 2), the axial load is immediately removed (unloading rate = 10 μm/s) and then the confining pressure is incrementally reduced again. Once a yield curve has been probed over a range of confining pressure increments, the pressure is increased again beyond P* I so that the sample sits on another new yield curve ( Figure 2a, step 3, P* II ). The procedure of incremental confining pressure reduction and axial loading is then repeated for this subsequent yield curve ( Figure 2a, step 4). This loading history therefore maps out the family of yield curves associated with the accumulation of purely isotropic volumetric strain (ε 1 = ε 2 = ε 3 ). All the movement in stress space is inside a given yield envelope in the elastic regime; the only time that this is not the case is at the onset of yield itself. Typically, in order for yield to be identified during a given stress-probing increment, the sandstones are deformed beyond yield to an axial strain of less than 0.1%. This is comparable to the axial strain required to identify yield in bassanite in the study of Bedford et al. (2018; see their Figure  6c). Bedford et al. (2018) showed that any inelastic damage produced during this stress-probing is minor and has little effect on the yield curve being mapped (see Figure S3 in their supporting information).

Deviatoric Overconsolidation
The second type of overconsolidation test that can be performed in a traditional triaxial configuration is deviatoric overconsolidation (axisymmetric compression). Unlike the hydrostatic overconsolidation tests, this loading path has a component of inelastic shear strain as well as volumetric strain. The loading history is summarized in Figure 2b. Initially, a sample is loaded to a given pressure below P* before being axially loaded to induce some inelastic deviatoric strain ( Figure 2b, step 1). After a given increment of deviatoric loading, the load is removed to move back into the elastic regime and the confining pressure is incrementally increased or decreased so that the yield curve can be probed ( Figure 2b, step 2). Following this, the confining pressure is returned to the value used for the initial axial loading and the sample is further deformed along the same loading path (Figure 2b, step 3). The loading is then stopped again after another increment of inelastic deviatoric strain has been accumulated and the probing procedure is repeated (Figure 2b, step 4).

Microstructural Analysis
Upon removal from the deformation apparatus, samples were prepared for microstructural analysis. Samples were vacuum impregnated with epoxy resin while still in their polyvinyl chloride jackets to Samples are hydrostatically loaded beyond P* so that inelastic volumetric strain is accumulated and the sample sits on a new yield curve (step 1). The pressure is then incrementally reduced to move back inside this new curve and the sample subjected to elastic axial loading between each pressure reduction step (step 2). The axial load is applied until the deviation from linear elastic loading is observed (inset) upon which it is immediately removed. The figure shows some schematic probing increments with yield points (filled circles); however, for neatness, not all probing loading paths are shown (empty circles). Once the yield curve is fully probed, the sample is further overconsolidated so that more inelastic volumetric strain is accumulated (step 3). The subsequent yield curve is then probed in the same way to determine if there is any shape evolution (step 4). (b) A schematic summary of the deviatoric overconsolidation experimental technique. The sample is hydrostatically loaded to a given pressure and then subject to axial loading to accumulate an increment of inelastic deviatoric strain (step 1). The load is then removed and the sample subject to axial probing at different confining pressures (step 2). Once the new yield curve has been probed, the sample is further deviatorically overconsolidated along the same loading path so that more inelastic shear strain is accumulated (step 3). The subsequent yield curve is then probed in the same to monitor any shape evolution (step 4). ensure no material was lost. After impregnation, samples were cut and this surface was then reimpregnated to ensure resin fills as much of the pore space as possible. The samples were then polished ready for imaging in the scanning electron microscope (SEM). Backscatter electron (BSE) images were acquired using a Philips XL30 tungsten filament SEM using an accelerating voltage of 20 kV and a beam current of 3 nA at a working distance of 13 mm.

Initial Yield of the Starting Material
A list of traditional experiments performed to constrain the yield behavior of the starting material using multiple samples are presented in Table 2. Data on porosity evolution with increasing effective pressure are shown in Figure 3 with a clear deflection in the loading curves marking P* and the onset of inelastic compaction. There is some sample variability in the pressure that P* occurs; however, typical values are 37-48 MPa for Boise sandstone and 54-59 MPa for Idaho Gray sandstone (range determined from five samples per rock type). It can be seen that, from tracking the porosity evolution during pressure reduction, most of the deformation beyond P* was inelastic as the porosity does not return to the initial starting porosity once the sample has been unloaded.
A compilation of mechanical triaxial data for the starting materials, collected in the traditional method using multiple samples, is presented in Figure 4. In this study, we adopt the convention that compressive stresses and compactive strains are positive (i.e., axial shortening and porosity reduction). Figures 4a and 4b show differential stress versus axial strain data for multiple samples of Boise and Idaho Gray sandstones loaded over a range of effective pressures. Yield is marked by the deviation from quasi-linear elastic loading. Beyond yield, the majority of samples display strain hardening; only samples loaded under the lowest effective pressures (5 and 10 MPa) display postyield weakening. During strain hardening, some of the Idaho Gray samples experience transient stress drops (Figure 4b), which have been shown previously to be associated with compaction localization . However, upon inspection at the end of the experiment, no evidence for localization was observed. This might be because the relatively high strains to which the samples were deformed destroys any evidence of localization. It may be observed better if the deformation Note. This was done in the conventional sense using multiple samples.

10.1029/2018JB016719
Journal of Geophysical Research: Solid Earth were stopped and the sample examined immediately after a stress drop, however, this would require many extra tests and is beyond the scope of the current study. It is interesting to note that similar stress drops are not observed for Boise sandstone. It is not immediately apparent why Idaho Gray sandstone show stress drops potentially linked to localization, whereas Boise sandstone with very similar petrographic properties (Table 1) does not. Figures 4c and 4d show the effective mean stress versus porosity reduction during axial loading for Boise and Idaho Gray sandstones, respectively. In a conventional axisymmetric triaxial test, where the confining pressure is being held constant during axial loading, as the axial stress is increased by an increment Δσ 1 , the effective mean stress will increase by an increment of Δσ 1 /3. Hence, the porosity evolution with effective mean stress during axial loading can be compared to the porosity evolution from hydrostatic loading (Figures 4c and 4d). The data show that all samples experience shear-enhanced compaction once yield has been achieved, apart from the sample of Boise sandstone loaded from an effective pressure of 5 MPa, which experienced dilation. Figure 4e shows yield curves, plotted from the triaxial data, for both sandstones. The curves have a relatively shallow geometry compared to previous studies on sandstone (e.g., Wong et al., 1997) and, as with previous investigations employing this methodology, there is some data scatter of yield points around each curve. It is interesting that only the Boise sample loaded from an effective pressure of 5 MPa experienced dilation, as there are other yield points on the low-pressure side of the yield curve for both rock types. This may be highlighting general shortcomings in the assumptions that are typically made about the mode of deformation expected at different points along the yield curve.

Hydrostatic Overconsolidation Results (Single Samples)
A list of single sample overconsolidation tests can be found in Table 3. Yield curves from the hydrostatic overconsolidation tests are presented in Figure 5, and for comparison, the initial yield data presented in Figure 4e have been included (empty circles). The data show that the curves of both sandstones are very similar with a broadly elliptical shape, albeit a shallow one. The compactive side of the yield curves are also partly comprised of near-vertical limb in the vicinity of P*, similar to those reported by Bedford et al. (2018) and Tembe et al. (2007). Each curve in the family is labeled with the porosity at its respective P*. However, these curves are not strictly isoporosity contours as there is elastic unloading of the sample as the pressure is incrementally reduced during the stress-probing procedure. The total amount of relaxation that occurs as the sample is completely unloaded from P* is typically about a 2% porosity increase.
To test whether there is evolution in the family of yield curves, each curve can be normalized by dividing P and Q by the respective P* value ( Figure 6). This normalization technique has been used previously to compare yield curves from different rock types (e.g. Cuss et al., 2003b, Wong et al., 1997 as well as evolution with inelastic strain (e.g., Bedford et al., 2018). The normalized data for both sandstones reveal that there is only minor evolution of the yield curves in response to hydrostatic overconsolidation (i.e., increasing volumetric strain). The yield curves appear to expand slightly on the compactive side. Normalization further highlights the flattened shape of the yield curves. The peaks of the curves have a Q/P* ratio of 0.22-0.3. This is substantially lower than values that have been previously reported for other sandstones that typically have Q/P* ratios in the range of 0.5-0.7 (e.g., Wong et al., 1997). For comparison, data from the multisample tests to determine the initial yield curve of the starting materials (Figure 4e) have been included in Figures 5 and  6. These data show that the probing technique produces much higher resolution yield curves than those mapped by the traditional method using multiple samples, which have considerable data scatter (e.g., Figure 4e).
The validity of the stress-probing methodology used in this study was tested on Idaho Gray sandstone by deforming a fresh sample (IG16 and IG17 in Table 3) directly to the maximum amount of overconsolidation and comparing the results obtained in Figure 5 (Figure S2). The data show that there is small difference in magnitude between the different samples but the overall shape evolution is similar to that observed in

Deviatoric Overconsolidation Results (Single Samples)
Yield curves from the deviatoric overconsolidation tests are presented in Figures 7c and 7d with the dashed arrow marking the loading trajectory along which the sample was deformed. For comparison, these data are  Note. In total, 503 yield points were identified in the tests.

10.1029/2018JB016719
Journal of Geophysical Research: Solid Earth presented alongside yield curve evolution mapped in the traditional sense where individual yield points are obtained from different samples (e.g., Baud et al., 2006). As the samples from our multisample tests (section 3.1) were deformed considerably beyond their inelastic yield, we can use the volumetric strain data presented in Figure 4 to plot yield curve evolution in this traditional manner (Figures 7a and 7b) such that it can be compared to yield curves obtained from the single sample stress-probing methodology (Figures 7c  and 7d). The single-sample yield curves were mapped along the deviatoric loading path after the same amounts of inelastic volumetric strain (ε v ) as the yield curves plotted using multiple samples (1%, 3%, and 5%) so that direct comparisons between the methodologies can be made. The single samples were then deformed to greater volumetric strain along the same stress path in order to analyze further the effects on yield curve evolution.
A significant benefit of the single sample tests is that they allow for an almost entire reconstruction of the yield curve after a given amount of inelastic strain. This enables the observation that there is a significant expansion of the compactive side of the yield curve, with P* migrating to higher effective mean stresses with continued deformation (Figures 7c and 7d). It should be noted that P* was not always probed for all deviatoric strain increments as it can become difficult to identify unless the pressure is increased significantly beyond P*. Therefore, to ensure that no inelastic strain was accumulated during the stress-probing procedure, the confining pressure was conservatively increased to ensure that P* was not greatly exceeded. Despite this, the evolution of P* toward greater effective means stresses can still be inferred from the rest of the yield curve. A similar phenomenon was seen by Bedford et al. (2018) for porous bassanite (see their Figure 11). The same evolution, of an expanding P*, cannot be observed in the multisample tests (Figures 7a and 7b) as in these experiments, the data are only gathered along the stress trajectory of a given test. The reconstruction of an entire yield curve using this methodology is achieved by using the data from all the different samples deformed along different loading paths (dashed arrows in Figures 7a and 7b). The yield curve evolution will likely be different along all these stress paths; however, this cannot be constrained from an individual test as it is not known how the rest of the yield curve is evolving away from the given loading path and thus the complete evolution, such as an expanding compactive regime, will remain unconstrained. Figure 6. Normalized yield curves from the hydrostatic overconsolidation tests for (a) Boise and (b) Idaho Gray sandstones. Each curve is divided by its respective P* value. There is a slight expansion of the compactive side of the yield curves with increasing volumetric strain. Initial yield data for the starting materials that have also been normalized are included to highlight the difference in resolution between the different methodologies for mapping yield curves.
Only the experiments where loading was initiated from the same starting pressure (labeled in Figures 7a and  7b) can be directly compared to the yield curve evolution mapped from the single sample experiments (Figures 7c and 7d).
Another observation from the deviatoric overconsolidation tests is the appearance of a plateau at the peak of the yield curve, particularly at high strains (Figures 7c and 7d). This can also be seen in an additional deviatoric overconsolidation data set where the deviatoric loading was initiated at pressures greater than the P* values of the starting materials ( Figure S2 and Text S1). The emergence of a plateau at higher strains (ε v = 5-10%), as well as the significant expansion of the yield curve at these strains, shows that compaction is getting more difficult and the transition between the predominantly dilatant and the predominantly compactive regimes, which is typically considered to be the CSL (Figure 1b), develops into a broad zone rather than a line. The lack of P* data means that the deviatoric overconsolidation yield curves cannot be normalized in the same way as the hydrostatic overconsolidation curves. However, yield curve evolution in response to deviatoric strain can be analyzed by comparing with the hydrostatic overconsolidation tests as described in the following section.

Hydrostatic Versus Deviatoric Strain
As the hydrostatic overconsolidation tests only involve purely inelastic volumetric strain, any differences in the yield curve evolution during deviatoric overconsolidation is likely a result of inelastic shear strain as this loading path will impart both volumetric and shear strain on the sample. Comparisons can be made between yield curves from the different test type where samples have experienced a similar amount of volumetric strain (Figure 8). The data for both Boise and Idaho Gray sandstones show that there is a large expansion of the yield curve when samples are compacted along a deviatoric stress path. The peak stress of the deviatoric yield curves is approximately 3 times that of the hydrostatically compacted samples, even though the amount volumetric strain the samples have experienced is comparable. This highlights that inelastic shear strain has a substantial effect on the relative size of sandstone yield curves.  Figure 4. Ellipses are fitted to show how yield curves are mapped using the traditional multisample methodology. (c, d) Yield curves from the single sample deviatoric overconsolidation tests. The stress path is shown by the dashed arrow, and the yield curves were probed after the same amounts of inelastic volumetric strain as presented in panels a and b. The samples were then deformed further to accumulate more inelastic volumetric strain (ε v ): 6.5% for Boise sandstone and 7.2% for Idaho Gray sandstone.

10.1029/2018JB016719
Journal of Geophysical Research: Solid Earth

Microstructures
Undeformed samples of both sandstones show that the majority of quartz and feldspar grains have relatively few microcracks and fractures (Figures 9a and 9b). Postmortem examination of single samples used for the stress-probing experiments shows deformation occurred relatively homogeneously along both the hydrostatic and deviatoric loading paths, predominantly by cataclastic flow (Figures 9c-9f), with a lack of any obvious sample-scale structure or localization. Although the final porosities of the hydrostatically and deviatorically deformed samples used for microstructural analysis are different, meaning that a direct comparison cannot be made, the overall structure in both sets of images is broadly similar with distributed deformation being evident. It is clear that the greatest mechanical difference between hydrostatically and deviatorically compacted samples occurs at the highest strains. Therefore, it might be expected that any anisotropy responsible for the difference would be most evident in the microstructures from these high strain experiments (Figures 9c-9f). The high strains potentially could have eradicated earlier evidence of localization or any significant structure, but this still does not answer the question of what is responsible for the significant difference in the size of the yield curves. Tests could be conducted and stopped at lower strains to observe any development of the microstructure, but these would be related to smaller differences in the mechanics, and contrasts in the microstructures might be more difficult to distinguish.
From the microstructures observed, it is clear that the deformation is dominated by intragranular cracking that leads to grain crushing, which has been shown previously to be the dominant hydrostatic and triaxial compaction mechanism at the onset of yield in siliciclastic rocks . The deviatorically overconsolidated samples were compacted to a lower porosity and hence show more intense deformation (Figures 9e and 9f). Cataclasis and grain crushing leads to the formation of a poorly sorted fine-grained deformation product (in the order of a few microns), which anastomoses around large intact relict grains. It is likely that a large amount of the postyield deformation can be accommodated by grain rearrangement in this finer-grained material after it has formed, which leads to the preservation of large relict grains.  (c, d) Samples that have been hydrostatically overconsolidated with evidence for grain crushing (porosity = 16.8% and 17.5%, respectively). (e, f) Deviatorically overconsolidated samples with fine-grained deformation product anastomosing around relict large grains (porosity = 11.5% and 10.0%, respectively, samples from data presented in Figure S2). The maximum principal stress (σ 1 ) is orientated vertical with respect to the images.

Journal of Geophysical Research: Solid Earth
It has been documented in previous studies on sandstone compaction that grain crushing is facilitated by Hertzian fractures that radiate from impinging grain contacts because of local stress concentrations (e.g., Menéndez et al., 1996;Fortin et al., 2007;. Similar Hertzian fractures are seen in these sandstones (Figures 10a and  10b) with intragranular cracks propagating across larger grains of quartz and feldspar, ultimately leading to grain crushing and pore collapse. Mica grains that were trapped between quartz and feldspar grains deform via grain kinking and cleavage cracking ( Figure 10c).
As there is little noticeable difference in the micromechanics between samples deformed along the different loading paths, we investigated if any microstructural anisotropy had developed during deviatoric loading by measuring the orientation of microfractures along horizontal transects across the BSE images. All BSE images were acquired with σ 1 vertical with respect to the image. The orientation data show that for a hydrostatically overconsolidated sample, the fractures are orientated fairly randomly (Figure 11a), which is consistent with Hertzian fractures radiating from grain contacts in an isotropic sample-scale stress field. However, for the deviatorically overconsolidated sample, there is a preferred orientation of fractures aligned at an angle of less than ±40°to σ 1 (Figure 11b) suggesting that a microstructural anisotropy develops during deviatoric compaction.

Discussion
Before the implications of the results are discussed, it is worth comparing the yield curve data acquired using multiple samples (Figures 4 and  7a and 7b) with those acquired from the overconsolidation tests on single samples (Figures 5, 7c and 7d, and S2). It is clear from looking at the results of the single sample tests that the resulting yield curves are of much higher resolution than those mapped in the traditional sense using multiple different samples. Although using enough samples appears to average out the effects of sample variability and allow for approximate ellipses to be fitted to the data (Figures 7a and 7b), in fact, this does not give a true representation of yield curve evolution along a particular stress path. This is because each individual sample is loaded from a different starting pressure and only the evolution of the yield points along that stress trajectory are observed. The complete yield curve evolution for a given stress path can actually be quite different, as observed in the single-sample data (Figures 7c and 7d) where there is a significant expansion of the compactive side and migration of P* to higher effective pressures. Therefore, only the single sample method will provide a complete picture of how the entire yield curve is evolving along a given stress path. The effects of sample variability are not completely negated by the use of a single sample, as the chosen sample could represent an outlier in the data and produce a curve where the magnitude would differ slightly from the norm for that rock type. Evidence for this was seen in the validation tests where the yield curve shape was relatively consistent between samples but the magnitude varied ( Figure S3 and Text S2). However, regardless of any slight differences in magnitude, the overall geometry of the yield curve mapped for a single sample is precisely constrained and any subsequent evolution with increasing inelastic strain can be determined relative to this curve.

Yield Curve Shape Evolution
A standout feature from the sandstones used in this study is that the yield curves have a very flat geometry in comparison to previous investigations. Wong et al. (1997) suggested that normalized yield curves Figure 10. (a, b) Quartz and feldspar grains deform by grain crushing during deviatoric loading, which is facilitated by the propagation of Hertzian fractures from grain contacts. The maximum principal stress (σ 1 ) is orientated vertical with respect to the images. (c) Mica grains deform by grain kinking and become pinched between quartz and feldspar grains.
should have a peak Q/P* ratio in the range of 0.5-0.7; however, normalized curves of the sandstones in this study have a peak Q/P* ratio of approximately 0.3 ( Figure 6). This low Q/P* ratio is seen for both the yield curves mapped using multiple samples (empty circles in Figure 6) and those from the singlesample hydrostatic overconsolidation tests, suggesting that the low aspect ratio of these yield curves is not a result of the stress-probing methodology. Both Boise and Idaho Gray sandstones have porosities that are higher than any of the sandstones reported by Wong et al. (1997) and are therefore probably more comparable with loose or poorly lithified sands. However, normalized data for loose sand (e.g., Crawford et al., 2004;Miura et al., 1984;Nguyen et al., 2014) also typically fall in between the bounds proposed by Wong et al. (1997). The similarities in yield curve shape between loose sand and lithified sandstone rule out the strength of grain contacts as the reason for the flatter yield curves of this study. Crawford et al. (2004) investigated the effects of sand textural characteristics on yield curve shape. They found that Ottawa Sand, which is consists of coarse (approximately 1 mm) well-rounded grains of quartz, had a relatively shallow yield curve and was the only sand they tested that when the data were normalized plotted below the lower bound of Wong et al. (1997), with a Q/P* ratio just below 0.4. Both Boise and Idaho Gray sandstones have a similar coarse grain size to Ottawa Sand, which could potentially explain the similar flatter yield curves; however, they differ texturally in that the grains are more angular. One major difference between the sandstones of this study and many investigations of loose sands is the modal composition. Ottawa sand and many of the other loose sands that have been studied are predominantly comprised of quartz, whereas both Boise and Idaho Gray sandstone are only about 50% quartz with a large amount of feldspar and other minerals (Table 1). It has been shown previously that mineralogy can affect the size of yield curves, but its effects on shape remain unclear. For example, Bleurswiller sandstone has a high abundance of minerals other than quartz (e.g., feldspar, mica, and oxides) and its compactive yield stresses are lower, by approximately a factor of 5, than Bentheim sandstone, which is similar porosity but comprised almost entirely of quartz (95%; Tembe et al., 2008). However, the normalized peak Q/P* ratios of these sandstones are broadly similar (approximately 0.45 for Bleurswiller and 0.53 for Bentheim), and it remains uncertain whether composition can also lead to a systematic variation in yield curve shape as well as size. Perhaps the flatter yield curves observed for Boise and Idaho Gray sandstones in this study are a combined result of both texture and mineralogy. Irrespective of the reason for the eccentric yield curves, the data suggest that caution should be taken when generalizing yield envelopes (e.g., Rutter & Glover, 2012;Wong et al., 1997;Zhang et al., 2000) and that the previously ascribed bounds for normalized yield data could be extended to incorporate yield curves with a flatter geometry, particularly for overconsolidated porous rocks.
As well as the yield curves having a flattened geometry, the data also reveal that there is evolution in shape with continued deformation. The hydrostatic overconsolidation tests show a slight expansion of the compactive side of the yield curve with increased volumetric strain (Figure 6). A similar yield curve evolution with increasing volumetric strain was also observed for the sandstones studied by Baud et al. (2006). However, Figure 11. Rose diagrams of microfracture orientation from (a) hydrostatically overconsolidated and (b) deviatorically overconsolidated samples of Idaho Gray sandstone. The hydrostatically deformed sample shows a distributed set of fracture orientations. The deviatorically deformed sample has fractures that are preferentially aligned with their long axis less than ±40°from σ 1 .

10.1029/2018JB016719
Journal of Geophysical Research: Solid Earth this is different to the evolution observed for porous bassanite in the study of Bedford et al. (2018) where the peak Q/P* ratio decreases with continued hydrostatic compaction. This is likely due to differences in the micromechanical response between bassanite and sandstone. Porous bassanite is made of platy grains, has a fracture-like pore network, and compacts predominantly via grain rearrangement and closure of these fracture-like pores, whereas sandstone is made of more equant grains and pores and compacts predominantly through grain crushing and pore collapse (e.g., Figure 10). Despite there being some yield curve shape evolution as a result of hydrostatic strain in this study, it is minor and the overall yield curve geometry is well constrained to a shallow ellipse with a steep limb on the compactive side. Bedford et al. (2018) described this near-vertical section of the yield curve as a region of "shear-insensitive compaction" and a similar region can be defined for Boise and Idaho Gray sandstones (Figure 12a).
Although yield curve evolution in response to hydrostatic compaction is relatively minor for both Boise and Idaho Gray sandstone, their response to deviatoric compaction is more significant. It is observed in the deviatoric overconsolidation tests that a plateau forms at the peak of the yield curve as more inelastic shear strain is accumulated (Figures 7c, 7d, and S2). The plateau is particularly prevalent once 5-10% volumetric strain has been achieved and also coincides with a significant expansion of the yield curve. This suggests that along the deviatoric loading path, there may be two regimes of yield curve evolution; at low to moderate strains (ε v = 1-5%), the yield curves retain their curvature and any evolution is moderate, whereas at high strains (ε v > 5%), there is the formation a plateau and significant expansion in yield curve size. Similar plateaus were also observed by Bedford et al. (2018) on porous bassanite. They termed this region "pressure insensitive deformation" (Figure 12b) and associate it with a greater contribution of dilational processes to the bulk deformation behavior of the sample, even if the sample is still experiencing bulk compaction. It is known already that during deformation of porous rock both compaction and dilation occur together (e.g., Brace, 1978). During compaction of siliciclastic rocks, grain crushing is facilitated by the propagation Figure 12. Interpretation of the yield curves mapped in the hydrostatic and deviatoric overconsolidation experiments. (a) The curves from the hydrostatic overconsolidation tests are partly comprised of a near-vertical limb on the high-pressure side, which defines a region of shear-insensitive compaction. (b) The curves from the deviatoric overconsolidation tests show the appearance of a subhorizontal plateau as the sample accumulates more inelastic shear strain. This marks a transitional region of pressure-insensitive deformation, which is characterized by simultaneous compaction and dilation processes.

10.1029/2018JB016719
Journal of Geophysical Research: Solid Earth of Hertzian fractures across individual grains (e.g., Figure 10), which involves the opening of voids and hence dilatancy on the grain scale. As a sample experiences further compaction and porosity is reduced, there becomes less space available for grain rearrangement, and therefore, dilatancy can become a more dominant process in the bulk behavior of the sample. The formation of a plateau could therefore arise when compactive and dilational mechanisms reach a similar order as a sample deforms along the deviatoric loading path (Figure 12b). An alternative explanation could be that there is some sort of memory effect related to the stress path. For instance, during deviatoric loading, the sample is deformed to a given axial strain and the differential stress value obtained could potentially exert a greater control on any subsequent deformation than the effective mean stress. Therefore, any further deformation would only occur once this value is exceeded, even under different effective mean stresses, and this would produce a plateau-like feature on the yield curve (analogous to the Kaiser effect (e.g., Lavrov, 2003)). More detailed investigation is required to elucidate the micromechanisms responsible for the formation of a plateau. However, regardless of the explanation for its appearance, the critical state model, where the yield ellipse is divided into a dilatant low-pressure side and a compactive high-pressure side separated by the CSL (e.g., Figure 1), is an oversimplification. Rather, the CSL develops into a critical state zone with the mode of deformation likely occurring as a spectrum from dilation to compaction around the yield curve rather than separated into discrete zones. Therefore, any attempt to make generalizations about the CSL (e.g., Rutter & Glover, 2012) should be made with caution. Wood and Maeda (2008) have also suggested that a similar critical state zone can develop during sand deformation as the grading and particle size of the sand change in response to particle crushing.

The Effect of Stress Path on Yield Curve Evolution
As well as changes in yield curve shape that occur during inelastic deformation, it is important to consider how different stress paths affect the relative magnitude of the yield curves. Figure 8 highlights that deformation along the deviatoric loading path, which involves the accumulation of both inelastic shear and inelastic volumetric strain, produces a yield curve with a peak that is considerably higher than a sample that has been hydrostatically compacted (i.e., purely inelastic volumetric strain) to a similar amount of inelastic volumetric strain. This difference must be caused by the effects of inelastic shear strain along the deviatoric overconsolidation pathway and is particularly pronounced at higher strains. Similar behavior was documented by Bedford et al. (2018), which they attributed to the development of a heterogeneous microstructure and conjugate shear bands. Although there is no obvious structure that develops in the Boise and Idaho Gray sandstones, it is likely that the reason for the difference in yield curve magnitude is microstructural anisotropy that develops along the deviatoric loading path. For other materials such as unconsolidated soft clays, microstructural anisotropy has been shown previously to exert significant influence on the size and inclination of yield curves (e.g., Sultan et al., 2010;Wheeler et al., 2003). For granular materials it has been shown that anisotropy in intragranular fractures, with a preferred orientation aligned parallel to the maximum principal stress (σ 1 ), develops during loading under a deviatoric load (e.g., Kanaya & Hirth, 2017, Wu et al., 2000. Figure 11b shows that there is a preferential alignment of microfractures with respect to σ 1 in the deviatorically overconsolidated samples of our study. Recent works have highlighted the importance of pore geometry and orientation on rock strength (e.g., Bubeck et al., 2017, Griffiths et al., 2017, with samples containing elongate pores (i.e., cracks and fractures) being strongest when compression is applied parallel to the long axis of the pores, compared to compression applied parallel to the short axis (Bubeck et al., 2017). Therefore, a set of preferentially orientated microfractures subparallel to σ 1 that develop during deviatoric loading might explain the apparent increase in strength and increase in yield curve magnitude relative to hydrostatically overconsolidated samples. If these tests were performed Figure 13. Schematic summary of the effect of loading path on yield curve evolution. Curves of the same color have been compacted to the same porosity (Φ) by (a) hydrostatic compaction and (b) deviatoric compaction. Deviatoric loading produces yield curves that are greater in magnitude than those formed by hydrostatic loading. A plateau also forms at the peak of the curves in response to deviatoric compaction. in a true triaxial apparatus, it may be possible to test the yield curve evolution in different directions across a deviatorically compacted sample to see if a mechanical anisotropy has in fact been induced by the formation of preferentially orientated microfractures. Future studies could also investigate the extent of anisotropy by measuring P wave velocity evolution along different stress paths. The effects of different loading paths on yield curve magnitude are summarized in Figure 13.
The results of this study emphasize that a complete understanding of the petrophysical controls on porous rock deformation and their influence on micromechanics is still lacking. Future studies should aim to systematically investigate how different material parameters (e.g., grain size/grain size distribution, mineralogy, and grain shape) affect deformation behavior. It seems clear that a universal model for yielding behaviour, based on porosity and/or grain size (e.g., Figure 1), is an oversimplification and more effort needs to be directed toward factors that influence the evolution of yield curves with inelastic strain. Particularly for high strains, for example, those that might be occur locally around the wellbore region of a reservoir, as this is when significant evolution can occur (Figures 7c and 7d). Evolution of sandstone strength and microstructure under different loading conditions has important implications for reservoir production. For example, depletion-induced compaction can lead to changes in the stress state, which will in turn affect the deformation behavior of the reservoir. It is often assumed that reservoir compaction follows a uniaxial stress path; however, studies have shown that during production, the stress states can evolve and the rock is rarely under uniaxial conditions (e.g., Hettema et al., 2000, Rhett & Teufel, 1992.

Conclusions
A series of stress-probing experiments was performed on Boise and Idaho Gray sandstones to investigate the shape and evolution of their yield curves. The yield curves were found to have broadly elliptical geometry, albeit much flatter than typically associated with porous sandstone (e.g., Wong et al., 1997). This is likely a result of the textural and mineralogical characteristics of the sandstones in this study. The compactive side of the yield curves is partly comprised of a near-vertical limb in the vicinity of P* in agreement with previous studies (Bedford et al., 2018;Tembe et al., 2007).
The evolution of the yield curves with work hardening is dependent on the nature of inelastic strain. Hydrostatically overconsolidated samples, experiencing purely volumetric strain, showed only minor yield curve shape evolution with a slight expansion of the compactive side of the curve. However, samples that were deviatorically overconsolidated, with a component of inelastic shear strain as well as inelastic volumetric strain, displayed a much more significant yield curve evolution. A plateau forms at the peak of the curve with continued deviatoric loading, likely as a result of dilational deformation becoming more important in the bulk behavior of the sample. Deviatoric deformation also produces yield curves that are greater in magnitude than similar porosity samples that have been hydrostatically compacted. This is likely a result of microstructural anisotropy that develops along the deviatoric loading path because of the presence of microfractures preferentially orientated with their long axis parallel to the maximum principal stress. Future work should consider the effects of inelastic shear strain on yield curve evolution and rock strength. It is also important to better constrain the microphysical controls on yield curve shape as the results of this study have shown that sandstones can have curves with a much more eccentric shape than previously reported.