High‐Resolution Mapping of Yield Curve Shape and Evolution for Porous Rock: The Effect of Inelastic Compaction on Porous Bassanite

Porous rock deformation has important implications for fluid flow in a range of crustal settings as compaction can increase fluid pressure and alter permeability. The onset of inelastic strain for porous materials is typically defined by a yield curve plotted in differential stress (Q) versus effective mean stress (P) space. Empirical studies have shown that these curves are broadly elliptical in shape. Here conventional triaxial experiments are first performed to document (a) the yield curve of porous bassanite (porosity ≈ 27–28%), a material formed from the dehydration of gypsum, and (b) the postyield behavior, assuming that P and Q track along the yield surface as inelastic deformation accumulates. The data reveal that after initial yield, the yield surface cannot be perfectly elliptical and must evolve significantly as inelastic strain is accumulated. To investigate this further, a novel stress‐probing methodology is developed to map precisely the yield curve shape and subsequent evolution for a single sample. These measurements confirm that the high‐pressure side of the curve is partly composed of a near‐vertical limb. Yield curve evolution is shown to be dependent on the nature of the loading path. Bassanite compacted under differential stress develops a heterogeneous microstructure and has a yield curve with a peak that is almost double that of an equal porosity sample that has been compacted hydrostatically. The dramatic effect of different loading histories on the strength of porous bassanite highlights the importance of understanding the associated microstructural controls on the nature of inelastic deformation in porous rock.


Introduction
Fluids are abundant in the Earth's crust and are of utmost geological importance for a variety of natural and industrial processes. These range from the transport of chemical components in metamorphic and hydrothermal, possibly mineralizing, systems (e.g., Bedford et al., 2017;Plümper et al., 2016) to the expulsion or extraction of hydrocarbons from a reservoir. Any fluid phase in the Earth's crust must be stored in some form of available pore space. The nature of this pore space and its evolution in response to tectonic stresses will in turn determine the ability of a rock to store and transport fluids. The response can be purely mechanical where the rock will either compact or dilate depending on the stress conditions, or porosity can change as a result of chemical processes, such as pressure solution and cementation. This study will focus on the initial time-independent mechanical response of porous rock to macroscopic stresses; an understanding of which can provide important insights into faulting and fluid flow in porous formations (e.g., Aydin & Johnson, 1978), as well as industrial problems such as borehole stability (e.g., Dresen et al., 2010) and reservoir compaction (e.g., Fisher et al., 1999;Nagel, 2001).
As a rock is put under load it will have an initial quasi-linear elastic response (Figure 1a), meaning that any deformation is recoverable if the load is removed. However, if the stress is increased further, and the yield stress is reached, then the rock will move out of the elastic regime and accumulate permanent inelastic strain. The stress at which the onset of inelastic strain occurs is called the yield point and can be identified by the deviation from linear-elastic loading (Figure 1a). Materials that deform by brittle processes (such as compaction and dilation) are pressure-sensitive, meaning that the stress at which yield occurs varies with pressure. A yield curve can therefore be defined in stress-space to define the locus of points at which yield occurs. For porous materials the yield curve is typically plotted in P-Q space, where P is the effective mean stress (P ¼ P À P f where P ¼ σ 1 þ σ 2 þ σ 3 3 is the mean stress and P f is the pore fluid pressure) and Q is 1. Yield in porous materials can also be achieved under a purely hydrostatic stress state (i.e., no differential stress). The yield curve therefore intersects the P axis at the hydrostatic yield point, which is commonly referred to as P* and marks the onset of grain crushing and pore collapse (Zhang, Wong, & Davis, 1990b). Much of our understanding on the yield of porous rock is based on principles from critical state soil mechanics (e.g., Schofield & Wroth, 1968;Wood, 1991). Soil mechanics has shown that yield curves are typically elliptical in shape ( Figure 1b) and many empirical data on a variety of porous rocks such as sandstone (Baud et al., , 2006Baud, Zhu, & Wong, 2000;Cuss et al., 2003;Louis et al., 2009;Wong et al., 1997), limestone Baud, Schubnel, & Wong, 2000;Cilona et al., 2014;Vajdova et al., 2004), volcaniclastics (Zhu et al., 2011), and dehydrated serpentinite (Rutter et al., 2009) have been in general agreement with broadly elliptical-shaped curves. However, rocks are not soils, so there is no fundamental reason that yield curves should be elliptical in shape as in soil mechanics. In fact, studies on sandstone Fortin et al., 2005) and porous volcanic rock (Heap et al., 2015) have reported that both the low-and high-pressure sides of the yield curve can be linear, with positive and negative slopes, respectively. Alternative models, based on granular mechanics (Guéguen & Fortin, 2013) or dual porosity (Zhu et al., 2010), have been used to describe linear yield curves. However, it remains unclear what the underlying controls are on yield curve shape for porous rock. The initial response is elastic with a linear stress-strain relationship. At a certain stress the rock will begin to yield as it starts to accumulate inelastic strain. This can be identified by the deviation from linear elastic loading. (b) Yield points can be identified at a range of effective pressures to define a yield curve. The curve is typically thought to be elliptical in shape with the low-pressure side being associated with localized dilatant behavior and the high-pressure side with distributed compaction. The point where the curve intersects the effective mean stress axis is referred to as P*. (c) A family of yield curves can be defined for different porosities. The peaks of these curves define the critical state line which separates the regions of dilation from compaction. (d) A 3D representation of the yield envelope where the third axis is porosity (Φ) (assuming constant grain size). The ellipses from (c) space out along the porosity axis and their respective P* points can be joined to form the normal consolidation line. The trajectory L-M, along the outside of the yield envelope, is a schematic representation of the P-Q-Φ path of the postyield experiments in Figure 4.

10.1002/2017JB015250
Regardless of the shape of the yield curve if the stress state falls inside the curve, the rock will be in the elastic regime and any deformation will be recoverable. If the stress conditions are outside or on the curve, then the rock will experience some form of inelastic deformation. The mode of inelastic deformation that a rock experiences at a particular point on the yield curve can be, in its simplest form, separated into two regimes (Figure 1b). At low pressures, the mode of deformation is dilatancy, which is typically associated with localized deformation and faulting (e.g., Menéndez et al., 1996). At high pressures a rock will undergo shear-enhanced compaction, which is associated with distributed cataclastic flow (Curran & Carroll, 1979). This transition from localized dilatancy to distributed compaction is often referred to as the low-temperature brittle-ductile transition in porous rock (e.g., Rutter & Hadizadeh, 1991;. As a rock accumulates inelastic strain and the porosity evolves, the yield curve will expand or contract depending on the nature of deformation (Figure 1c). A family of yield curves can therefore be defined for rocks of different porosity (and/or grain size) or for a single rock undergoing compaction. The crest of each yield curve, which separates dilation from compaction, can be joined by a line called the critical state line (CSL). This can also be visualized in 3D where the continuum of yield curves forms a yield surface (Figure 1d). In this diagram the third axis represents porosity (Φ) and the family of ellipses space out along the porosity axis. The respective P* values of the individual yield curves form another line called the normal consolidation line (NCL) (e.g., Rutter & Brodie, 1995), also referred to as the hydrostat.
The concept of a family of curves, where the shape of the curve remains the same with only the size changing according to porosity and grain size, has been applied widely to unconsolidated sediments. A series of complex probing tests on both clay (Graham et al., 1983) and sand (Miura et al., 1984;Tatsuoka & Ishihara, 1974) have revealed the shapes of these families and found that a unique curve shape can be defined for the different materials. However, the applicability of this concept to porous rock remains unclear and experimental data on curve evolution are sparse as the high pressures required to reach P* for porous rock make much of the postyield behavior unattainable in many deformation apparatus. Baud et al. (2006) tested yield curve evolution for a range of porous sandstones and found evidence that the shape of the compactive side of the curve evolves with inelastic volumetric strain. However, the data of Baud et al. (2006) are subject to a large amount of data scatter as a result of sample variability making it difficult to determine reliably the exact shape of yield curves in order to monitor their evolution. Alternatively, other studies have suggested that the yield curve for porous rock can be constrained to a shape that is independent of rock type (Rutter & Glover, 2012;Wong et al., 1997). These contrasting viewpoints highlight the current uncertainty in our understanding of the specific shape and evolution of yield curves for porous rock.
In this study the compaction behavior of dehydrated gypsum is tested to see if it fits into a relatively simple critical-state model with elliptical yield curves or whether another interpretation is needed. The data are important for understanding the compaction of rocks undergoing solid volume reductions during metamorphic reactions, but also have wider applications to the mechanics of porous rocks in general. Like the majority of devolatilization reactions, the dehydration of gypsum is associated with a large solid-volume reduction, which generates porosity. It is important to understand the porosity evolution in a dehydrating system because if fluids are not able to drain efficiently, then fluid overpressures may develop, which could generate earthquakes or slow slip in subduction zone settings (e.g., Gao & Wang, 2017;Hacker et al., 2003). A soil mechanics approach has previously been used to investigate the deformation of porosity created as a result of metamorphic devolatilization reactions (Rutter et al., 2009;Rutter & Brodie, 1995); however, a detailed investigation of the compaction behavior of dehydrated gypsum has been lacking. Gypsum has been widely used to study processes associated with dehydration reactions as it provides an ideal analogue material (e.g., Bedford et al., 2017;Fusseis et al., 2012;Ko, Olgaard, & Wong, 1997;Leclère et al., 2016;Llana-Funez et al., 2012). The complete dehydration of gypsum, to form the product mineral bassanite, generates a porosity of 27-28%, which is comparable to many of the sandstones used in previous compaction studies (e.g., Cuss et al., 2003;Fortin et al., 2005). Understanding the compaction of this porosity is important for gaining a better knowledge of the interplay between fluid flow, deformation, and reaction during gypsum dehydration, which in turn gives insights into dehydration reactions in general.
After this introduction the results of this study are presented in two sections. In section 2 an initial set of conventional triaxial experiments are performed to test the yield and postyield behavior of porous bassanite to see whether the elliptical-shaped critical-sate model (Figures 1b-1d) is applicable to this material.

10.1002/2017JB015250
However, these experiments highlight discrepancies between the compaction behavior of porous bassanite and the simple elliptical yield curve shape of the critical-state model. Consequently, in section 3, a series of novel stress-probing measurements were developed to try to map precisely the yield curve to explore its true shape and subsequent evolution. These tests reveal that the yield curve evolves significantly with inelastic strain and that the evolution is dependent on the nature of the stress path leading to this strain, with samples deformed under a deviatoric stress having much larger yield curves than samples deformed purely by hydrostatic compaction. The implications of this new data and the wider applications for porous rock deformation in general are then discussed in section 4.

Testing the Elliptical Critical-State Model for Porous Bassanite
The type of experiments that are performed to test the elliptical critical-state model are conventional drained triaxial tests, where a suite of samples are loaded to different effective pressures before an axial load is applied and the stress-strain evolution is monitored (i.e., axisymmetric compression). In this setup the confining pressure, pore pressure, and axial load can all be varied and controlled independently from each other. The pore fluid pressure was held constant at 20 MPa for all tests in this study (including those in section 3). A more detailed description of the experimental apparatus can be found in the supporting information (Text S1 and Figure S1).
The data presented in the rest of this section outline how the critical state model was tested for bassanite. First P* was identified and the NCL was determined by subjecting samples to hydrostatic loading (section 2.2). The rest of the yield curve of the starting material was then mapped by axially loading samples at different pressures below P* (section 2.3). The postyield envelope was then investigated (section 2.4) by loading samples that were already sat on the yield envelope at pressures greater than P* (path L-M on Figure 1d). The implications of the data and the reasons why the mechanical behavior of porous bassanite does not agree with an elliptical critical-state model is discussed in section 2.5.

Sample Preparation
The starting material is natural alabaster gypsum from Volterra, Italy, which has been widely used in a range of studies investigating the dehydration and deformation of gypsum (e.g., Bedford et al., 2017;Fusseis et al., 2012;Ko et al., 1997;Leclère et al., 2016;Llana-Funez et al., 2012). Volterra gypsum is fine grained (50-200 μm), with an initial porosity of <1.5%, and is considered to be relatively isotropic; however, weak shape-preferred orientations of grains have been reported (Hildyard et al., 2011). Cores (50 mm length × 20 mm diameter) were drilled in the same orientation from a block of Volterra gypsum and were then precision-ground to square the ends to tolerance (±0.01 mm). Once the gypsum had been cored and squared, samples were placed in an oven where they were left to dehydrate to form bassanite. The oven temperature was set to 80°C, and reaction progress was monitored by measuring the weight loss associated with fluid expulsion ( Figure S2). Samples would typically take about 30 days to dehydrate completely to form bassanite, and they would be left in the oven until they were used for compaction tests to minimize the possibility of rehydration at ambient conditions. Once removed from the oven, the samples were left for a few minutes to cool and then the porosity of the dehydrated material was measured using He-pycnometry. All of the reaction was therefore complete prior to insertion into the deformation apparatus.

Identification of P* and Mapping of the NCL
In order to investigate the nature of the yield curve for porous bassanite it is important to establish first the point at which yield occurs under hydrostatic conditions (P*) and also the subsequent porosity reduction with increasing pressure after P* has been reached. Once a sample was inserted into the apparatus, the confining pressure (P c ) was increased to 23 MPa and the pore fluid pressure (P f ) to 20 MPa. Both parameters were increased incrementally to ensure that an effective pressure of 3 MPa was never exceeded so that no inelastic compaction occurred. Once these starting conditions were reached, the servo-controlled confining pressure pump was used to increase continuously the confining pressure at a rate of 0.01 MPa/s while monitoring the change in pore volume with the pore pressure volumometer (porosity proxy). The confining pressure pump is able to increase the pressure in the pressure vessel by approximately 20-25 MPa before it reaches the end of its stroke. In order to map the rest of the hydrostat the pressure was then increased incrementally with a Journal of Geophysical Research: Solid Earth 10.1002/2017JB015250 manual air-driven hydraulic pump, and between each increment the volumometer was allowed to equilibrate before the pressure was increased again. Figure 2 shows the results of hydrostatic loading on the porous bassanite samples. A clear deflection in the curve, corresponding to the yield point (P*), can be seen during the initial loading with the servo-controlled pump at an effective pressure of approximately 6 MPa. This value is much lower than typical P* values for sandstones with comparable porosity to bassanite. For example, Penrith sandstone (Φ = 28%) has a P* of 140-170 MPa (Cuss et al., 2003) and Bleurswiller sandstone (Φ = 23-25%) has P* values in the range of 135-195 MPa (Fortin et al., 2005;Tembe et al., 2008), highlighting the relative weakness of porous bassanite compared to other rock types. However, this means that much of the postyield behavior of bassanite is attainable and the majority of the hydrostat in Figure 2 does in fact represent postyield deformation. This is highlighted further by unloading of the sample to reveal that the majority of the compaction has been inelastic.

Determination of the Yield Curve for the Initial Material
As P* is such a low value for porous bassanite (approximately 6 MPa) it is difficult to map the initial yield curve as there is only a small pressure range below P* in which to work. Traditionally, the yield curve of porous rock is determined by axially loading samples at different pressures below P* (e.g., Wong et al., 1997). Figure 3a displays mechanical data from axial loading of three different bassanite samples at effective mean stresses of 2, 3, and 4 MPa at a constant displacement rate of 0.5 μm/s (strain rate ≈ 10 À5 s À1 ). All three samples displayed transient strain hardening after they had yielded before ultimately displaying weakening, which was also associated with small stress drops. The pore volume data ( Figure 3a) reveal that all three samples experienced a switch from compaction to dilation during loading; however, the switch to bulk dilation always occurred in the strain hardening regime before the peak stress was reached. The maximum compaction at these low effective stresses is <6% volumetric strain, which is quite small in comparison to data on similar porosity sandstones, such as Bentheim sandstone which experienced up to 16% volumetric strain before a switch to dilation (Baud et al., 2006). This highlights further the relative weakness of porous bassanite in comparison to other rock types. Upon removal from the deformation apparatus it was observed that all three samples exhibited sample-scale shear fractures, which likely formed as the deformation behavior transitioned from hardening to weakening and the stress drops occurred.
As the samples were loaded at such low effective pressures it is difficult to identify accurately where yield occurs. Figure 3b displays a zoom on the initial parts of the loading curves where slight deviations from quasi-linear loading can be seen. We interpret these deflections to be the yield point and mark the onset of shear-enhanced compaction (C*). When plotted on a P-Q diagram (Figure 3c), the points, along with P*, are found to be in agreement with a broadly elliptical shaped yield curve as predicted by critical-state soil mechanics, although it is difficult to constrain accurately the true shape with such few data points.

Mapping of the Postyield Surface
To gain a more comprehensive insight into the yield behavior of porous bassanite, tests were performed to map out the yield surface at pressures above P*. Although this does not provide any extra information on the initial yield state of bassanite, it will elucidate to the postyield deformation and help reveal the shape of the family of yield curves (Figures 1c and 1d) as porosity decreases. In some studies, it has been implicitly assumed that postyield deformation tracks along a unique yield surface illustrated in Figure 1d. However, the nature of the inelastic deformation may fundamentally change this yield surface that may only be applicable for the initial yield of porous materials. Consequently, it is important to understand if the size and shape Figure 2. Hydrostatic loading of porous bassanite. The inset in the figure is a zoom of the initial loading of bassanite samples using the servo-controlled pump to increase the effective mean stress at a continuous rate of 0.01 MPa/ s. A deflection can be seen at an effective mean stress of approximately 6 MPa marking the onset of inelastic deformation. The dashed lines represent the loading history that is lost as the sample is pressurized to an initial effective mean stress of 3 MPa; the solid lines are therefore extrapolated back assuming linear elasticity. The rest of the hydrostat was mapped by incrementally increasing the effective mean stress using a manual pump. These incremental increases in pressure are marked by the points on the curve, to just above a mean stress of 200 MPa. The porosity was also tracked during unloading to display the extent of inelastic deformation.

10.1002/2017JB015250
of the yield surface changes with inelastic deformation achieved along various loading paths. These tests were performed by loading bassanite samples to different effective pressures greater than P*, before being subjected to axial loading at a constant displacement rate 0.5 μm/s. The starting effective pressures in the tests are higher than the 6 MPa value of P* meaning that they are already sat on the NCL before the axial load is applied. Therefore, as a sample is subject to a differential stress, the trajectory in P-Q-Φ space should traverse along the compactive side of the yield surface toward the critical-state line (path L-M on Figure 1d). Thus, the shape of the yield surface should be revealed.
The data from the postyield tests are summarized in Figure 4. All samples show initial strain hardening (Figure 4a), which is to be expected when in the compactive regime. As the samples are already sat on the yield surface there are no deflections in the loading curves that would typically be seen if loading started from the elastic regime inside the yield surface. The majority of strain induced by axial loading is therefore inelastic. Samples that were loaded at the lowest effective pressures (10, 30, and 50 MPa) went through a transition from shear-enhanced compaction to dilation (Figure 4b). Upon removal from the rig it could be seen that the deformation in these samples had localized onto a sample-scale shear band. This is likely associated with small stress drops observed in the stress-axial strain data (Figure 4a). Samples that were deformed at higher effective pressures (70 and 90 MPa) only experienced strain hardening and compaction.

Implications of the Data for the Shape of the Yield Surface
The data gathered from the postyield tests ( Figure 4) should provide enough information to map comprehensively the compactive side of the yield surface. The position of the critical-state line should also be revealed by the "hooks" in the porosity data ( Figure 4b). These hooks (often referred to as C* 0 ) have previously been observed in other studies and represent the change from compaction to dilation (e.g., Baud, Schubnel, & Wong, 2000;Vajdova et al., 2004Vajdova et al., , 2012. However, when the data are plotted in P-Q-Φ space such as in Figure 1d, there are peculiarities in the shape of the yield surface. These peculiarities are most clearly highlighted when the yield surface is observed along the differential stress (Q) axis and the data are projected down onto the P-Φ plane. In this view, the CSL is expected to lay roughly half way between the NCL and the Φ axis for elliptical yield curves ( Figure 5a). However, with bassanite, we find that what we interpret as the CSL, from the transition from compaction to dilation, lies almost directly on top of the NCL (Figure 5b). This suggests that the compactive side of the yield envelope should be near-vertical. As this contradicts the simple elliptical yield cap critical-state model, a new interpretation of the data is required. Two possible explanations could result in the CSL overlying the NCL. First the yield curves for these compacted materials could have an unconventional shape (i.e., not elliptical) specific to bassanite, which would have its peak overlying P*. Alternatively, the yield surface might not be a simple static surface, and perhaps the shape evolves as a sample accumulates permanent inelastic strain as suggested by Baud et al. (2006). The data presented in Figure 4 are not of the resolution to determine the first possibility. To investigate further, a different experimental technique needs to be adopted to constrain the complete shape of the yield curve after different increments of inelastic strain to identify if there is any evolution.

Investigating the Geometry and Evolution of the Yield Surface With Compaction
In order to investigate the complete shape of individual yield curves we perform novel experiments to probe a given curve multiple times at different places. To do this, a sample is overconsolidated and then the load is reduced to move back into the elastic regime (i.e., inside the yield curve). The sample is then taken to different confining pressures inside the yield curve and axially loaded to yield at each pressure. Using the principle that the yield point can be identified by the deviation from linear elastic loading (Figure 1a), the load is increased until this deviation is observed. Yield is much more easily identified than in the tests to determine the initial yield curve ( Figure 3b) as can be seen in Figure 6c. As soon as this deviation is seen, the axial load is immediately removed from the sample (unloading displacement rate = 10 μm/s) to ensure no permanent inelastic strain is accumulated. The confining pressure is then changed to a different value and the sample is subjected to another axial loading increment to probe the yield curve at a different point in stress space. All the movement in stress space is therefore inside a given yield envelope in the elastic regime. The only time this is not the case is at the onset of yield itself. Once one yield curve has been probed, the sample is overconsolidated further so that the sample sits on another new curve. This new curve is then probed in the same way to investigate if inelastic strain causes the shape to evolve.
To study comprehensively the effects of inelastic strain on the shape of the yield curve, different types of stress paths used to produce overconsolidation need to be considered. In a conventional triaxial experimental setup two types of stress path leading to permanent strain can be achieved: isotropic loading, which produces purely volumetric strain, and deviatoric loading, which has a component of shear strain. We therefore perform tests to analyze how these two types of loading paths affect the evolution of the yield surface. The first type of test is termed hydrostatic overconsolidation where the effects of purely volumetric strain are examined. The second type of test is deviatoric overconsolidation, which has a component of shear strain as well as volumetric strain. The experimental procedure and results of each type of test are outlined below.
3.1. Hydrostatic Overconsolidation 3.1.1. Experimental Procedure Figure 6a shows a schematic summary of the loading path for hydrostatic overconsolidation tests. Initially, a sample was loaded hydrostatically beyond P* so that it would experience inelastic volumetric strain along the NCL and be sat on a new yield envelope at lower porosity. The confining pressure was then decreased in increments; between each increment the sample was axially loaded to probe the shape of the new yield curve. In a triaxial deformation apparatus a typical axial loading path has a gradient of 3 when plotted in P-Q space. This is because as the axial load (σ 1 ) is increased the confining pressure (σ 2 and σ 3 ) is being held constant. Therefore, as the differential stress is increased by 1, the effective mean stress increases by 1 3 . As soon as the yield curve is reached (identified by the deviation from linear elastic loading) the axial load is quickly removed and then the confining pressure is reduced further. After the yield curve had been probed over a range of confining pressure increments the pressure was increased again beyond the new P* so the sample reaches another new yield curve. The steps of confining pressure reduction and axial loading were then repeated for this subsequent yield curve. This loading history will therefore map out the family of yield curves that form in response to purely volumetric compaction.
One concern with this technique is that the probing of the yield curve might result in small but significant amounts of permanent strain of the sample. Consequently, measurements were performed to test the probing technique and ensure that the assumed elastic loading history and incremental yielding do not impose any significant permanent deformation on the sample. In these tests the yield points were re-examined as the effective pressure was increased back toward P* inside the yield envelope following the first set of measurements. The two sets of yield points ( Figure S3) are found to be close to each other in stress space suggesting that the probing technique has a minimal effect on the mechanical properties of the sample. Small differences between the yield points could also be explained by hysteresis that might be expected as the pressure is decreased then increased again.

Results
The yield points that were collected during the repeated probing test are collated in Figure 6b. It can be seen that the use of a single sample has produced minimal data scatter helping to reveal in detail the shape of the yield curve. On first inspection of the yield curves in Figure 6b, the shape appears to be broadly elliptical. However, the compactive side of the curve is not entirely smooth and consists partly of a near-vertical limb. Each curve in the family is labeled with the porosity at its respective P*. However, these curves are not strictly iso-porosity contours as there is elastic unloading of the sample at lower effective mean stresses as can be seen in the porosity data in Figure 6c. The total amount of elastic relaxation that occurs as the sample is unloaded from P* is typically about a 2% porosity increase. To compare the shapes of different yield curves within the family, it is common to normalize the individual curves by dividing P and Q by their respective P* values (e.g., Cuss et al., 2003;Wong et al., 1997). If the shape of the yield curves obtained does not evolve with increased volumetric strain, it would be expected that they all would plot directly on top of each other when normalized. Figure 7 is a plot normalized curves from the hydrostatic overconsolidation tests. It shows that the peak of the curve decreases with respect to P* the more the sample is compacted. This demonstrates that there is not one unique yield curve shape that simply increases in size as porosity is reduced. It demonstrates that inelastic volumetric strain has a significant effect on the overall shape of the yield curve.

Deviatoric Overconsolidation 3.2.1. Experimental Procedure
Tests were also performed to investigate the effect of inelastic deviatoric strain on the shape of the yield curve. The loading history for this type of test is summarized in Figure 8a. These tests also involved initially Figure 5. (a) A schematic projection of a theoretical 3D yield envelope down on to the P-Φ plane. If the critical state line is projected onto this plane, assuming elliptical yield curves, it would be expected to lie half way between the normal consolidation line and the Φ axis. (b) Mechanical data from the bassanite deformation tests of this study. The axial loading data (red) sit directly on top of hydrostatic data (black). The hooks in the axial loading curves (C* 0 ) should mark the location of the critical state line as they represent the switch from compaction to dilation. The position of the hooks directly above the hydrostatic loading data suggests that the critical-state line would overly the normal consolidation line in contrast to that expected for elliptical yield curves.

10.1002/2017JB015250
loading the sample beyond P*, to an effective pressure of 50 MPa, so that the sample accumulated some inelastic volumetric strain and moved onto a new yield curve. The sample was then loaded axially (displacement rate = 0.5 μm/s) to induce deviatoric overconsolidation. After a given increment of deviatoric loading, the load was removed and the confining pressure was either increased or decreased incrementally so that the yield curve could be probed. Once the yield curve had been probed, the effective pressure was returned to 50 MPa, the value used for the initial axial loading, and the sample subjected to further deviatoric deformation along the same loading path. The loading was then stopped after another increment of inelastic strain had been accumulated and the yield curve was probed again. The sample was deformed along this loading path until 8% inelastic axial strain was achieved. The shape of the yield curve was tested after every 2% increment of inelastic axial strain along this loading path, where the axial strain is equal to the total axial strain minus the axial strain at the inflection point.

Results
The loading path and yield points from the deviatoric overconsolidation tests are displayed in Figure 8b. There are fewer yield points on the compactive side of the yield curves in these tests as P* is not known. The confining pressure was therefore conservatively increased when probing the compactive side of the yield curve to ensure P* was not exceeded hence preventing any further inelastic strain from being accumulated. As in Figure 6b, the curves are not iso-porosity contours as there is elastic unloading of the sample at low pressures. The labeled porosities in Figure 8b are from the highest pressure stress-probing increments for Figure 6. (a) A schematic summary of the hydrostatic overconsolidation experimental technique. Samples are hydrostatically loaded beyond P* so that inelastic volumetric strain is accumulated and the sample is sat on a new yield curve. 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. The axial load is applied until the deviation from linear elastic loading is observed 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. The subsequent yield curve is then probed in the same way to determine if there is any shape evolution. (b) Results from the probing technique for a single bassanite sample. The data points reveal that the family of yield curves are partly composed of the near-vertical limb on the high-pressure side. (c) Example of mechanical data from some probing increments. The data are labeled with the effective mean stress from which loading was initiated and the corresponding yield points are labeled in Figure 6b. The example data are all from the yield curve P* = 70 MPa. During each increment, the sample is loaded until yield can be identified (marked with a circle).

10.1002/2017JB015250
the respective curves. The obvious feature of these curves is that a near-horizontal "plateau" appears to form on top of the yield curve the more inelastic strain is accumulated. The gradient of this plateau decreases, and it flattens off the further along the loading path the sample is taken. Figure 9a shows a backscatter electron image of the fully dehydrated starting material prior to any deformation. When gypsum is dehydrated in an oven at atmospheric pressure, bassanite grains form with a platy morphology (Singh & Middendorf, 2007), which can make individual grains hard to identify. The platy grains appear to form clusters, which have a similar grain size to the gypsum protolith (50-200 μm). This could be as a result of crystallographic inheritance from the original gypsum, which has been suggested to occur in a previous study (Hildyard et al., 2011). However, adjacent clusters have a relatively random orientation with respect to each other (Figure 9a). The pore network is composed of quite elongate pores, similar to a crack-like porosity, which separate the individual bassanite grains and clusters.

Microstructural Analysis
Images of the hydrostatically and deviatorically overconsolidated samples are shown in Figures 9b and 9c, respectively. Both samples were compacted to a similar porosity (approximately 14.5%) so provide a useful Figure 7. Normalization of the yield curves from the hydrostatic overconsolidation tests by dividing each curve by its respective P* value. The peak of the curves decreases with respect to P* the more the sample is compacted and volumetric strain accumulated. Figure 8. (a) A schematic summary of the deviatoric overconsolidation experimental technique. Samples are hydrostatically loaded beyond P* and then subject to axial loading to accumulate an increment of inelastic shear strain. The load is then removed and the sample subject to axial probing at different confining pressures. 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. The subsequent yield curve is then probed in the same way to monitor any shape evolution. (b) Results from the deviatoric overconsolidation tests on a single bassanite sample. The loading path is marked by the black dashed arrow, and the yield curve was probed every 2% increment of inelastic axial strain. The data points reveal that the top of the yield curves starts to form a plateau the further along the deviatoric loading path the sample is taken.

10.1002/2017JB015250
representation of how different microstructures develop along different loading paths. The hydrostatically overconsolidated sample appears to have experienced fairly homogeneous deformation. Localized features are not easily identified. The nature of the bassanite grains makes it hard to identify if abundant grain crushing occurred during compaction. However, it is likely that the predominant deformation mechanisms are a combination of grain rearrangement and pore collapse. The deviatorically overconsolidated sample did localize deformation onto a sample-scale shear band as the strain approached 8%. A backscatter image away from the main shear band (Figure 9c) shows that there is also localized deformation in the body of the sample in the form of multiple conjugate sets of shear bands.

Discussion and Implications
It is clear from the initial tests conducted (section 2) that the critical-state line, in its traditional interpretation as the boundary between compactive and dilatant deformation, is not applicable. Rather, the data suggest that the transition occurs over a region and that within this region both compaction and dilation are occurring within the sample at the same time (e.g., Brace, 1978). This is highlighted further by the curvature of the strain hardening part of loading curves in Figures 3 and 4, and also that the transition to dilation occurs before the peak stress is reached (Figure 3). Consequently, the subsequent yield curve evolution will be affected by the trade-off between these competing processes. This could potentially make it possible to have a bulk transition from compaction to dilation (shown by the hooks in Figure 5) under stresses that would not be expected from the initial yield condition. Further evidence for this can be seen in the evolving yield curves of the hydrostatic and deviatoric overconsolidation tests, which will be discussed below.

Yield Curve Shape and Evolution
Data from both the hydrostatic and deviatoric overconsolidation tests have shown that yield curves are not perfectly elliptical in shape. The hydrostatic overconsolidation experiments have revealed that the compactive side of the curve is not smooth and partly composed of a steep limb. As this part of the curve is nearvertical it suggests that there is a region in P-Q space where yield is insensitive to differential stress. We term this region shear-insensitive compaction (Figure 10a). In this region an amount of differential stress needs to be overcome before significant shear-enhanced compaction will occur. Similar steep limbs have also been reported in a study where sandstone samples were axially deformed as pore volume was held constant so that the stress path would traverse along the compactive side of an individual yield curve (Tembe et al., 2007). An upsurge of acoustic emission activity was recorded by Tembe et al. (2007) as the samples moved from the steep limb to the more traditional curved part of the yield curve, which perhaps marks the onset of significant shear-enhanced compaction as suggested in this study.
The deviatoric overconsolidation tests have shown that the yield curve develops a plateau as the sample accumulates more inelastic strain. Similar plateaus have been hinted at using conventional triaxial tests on Darley Dale sandstone (Figure 4d from Baud et al., 2006). This near-horizontal portion of the yield curve is termed pressure-insensitive deformation (Figure 10b). In this region both compaction and dilational processes are expected to occur simultaneously and microstructural evidence for this will be discussed below in section 4.3. The concept of a critical-state line, separating a region of dilation from compaction (Figure 1c), could therefore instead be considered a "critical-state zone." The notion of a critical-state zone, rather than a line, has previously been suggested in the soil mechanics literature (Wood & Maeda, 2008). The formation of this zone is thought to occur as the grading and material properties of the soil change in response to particle crushing during deformation. Similar processes are in operation during porous rock deformation so it seems reasonable that a critical-state zone could also develop with continued pore collapse and grain crushing.
The shape of the yield curve evolves with inelastic strain along both the hydrostatic and deviatoric loading trajectories. Normalization of the hydrostatic overconsolidation data has shown that inelastic volumetric strain causes the yield curve shape to evolve (Figure 7). It has previously been suggested that the peak of a yield curve scales with P* (Cuss et al., 2003;Wong et al., 1997) and that this relationship can be used to define a unique curve independent of rock type (Rutter & Glover, 2012). However, the normalized data in Figure 7 highlight that there is variation in the ratio of Q/P* for the same sample as it accumulates inelastic strain. Although this is not a like-with-like comparison, as the previous studies were considering the initial yield condition, the data do suggest that caution should be taken when ascribing a universal yield curve to characterize general porous rock deformation. In particular, the changing aspect ratio observed in Figure 7 is likely a result of an evolving microstructure of the bassanite sample with continued deformation. There Figure 10. Interpretation of the yield curves mapped in the hydrostatic and deviatoric overconsolidation experiments. (a) The curves from the hydrostatic overconsolidation tests are partly composed of a near-vertical limb on the high-pressure side. This region is termed 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 is viewed as a transitional region of pressure-insensitive deformation, which is characterized by simultaneous compaction and dilation.
Journal of Geophysical Research: Solid Earth 10.1002/2017JB015250 will also undoubtedly be microstructural variations between different porous rock types making it likely that this will also cause variations in the aspect ratios of their respective yield curves.
The models of DiMaggio and Sandler (1971) and Carroll (1991) both use an evolving elliptical yield "cap" to model the deformation of porous sandstone with inelastic volumetric strain. In both cases the cap is used to describe the ductile compactive yield in conjunction with a different surface to model brittle (dilational) shear failure. This differs from the critical-state model (Schofield & Wroth, 1968) where the cap is part a single yield surface such as that in Figure 1b. Baud et al. (2006) have experimentally investigated the effects of inelastic volumetric strain on yield curve evolution for a range of sandstones. They found that the compactive side of the yield curve does evolve as a function of inelastic volumetric strain and that some sandstones appear to fit the DiMaggio and Sandler (1971) model, whereas others are better approximated by the Carroll (1991) model. However, the data collected by Baud et al. (2006) use many different samples from the same rock. This produces inevitable data scatter as a result of sample variability, which may hide subtleties in the true shape of the curve, such as a steep limb and a region of shear-insensitive compaction. Although the aspect ratios of the yield curves in the hydrostatic overconsolidation experiments of this study do change, neither model would fit particularly well with the data presented here as it does not appear necessary to use different surfaces to model the brittle and ductile parts of the yield curve. Rather, it would be more appropriate to use one curve to describe the regions of dilation and shear-enhanced compaction in conjunction with another curve to describe a region of shearinsensitive compaction.
The study of Baud et al. (2006) aside, there are a lack of previous experimental data on the effects of inelastic shear strain on yield curve evolution. Grueschow and Rudnicki (2005) incorporate the effects of inelastic shear strain into their constitutive model, which also considers the effects of inelastic volumetric strain. In this model, inelastic shear strain should cause the peak of the yield curve to expand relative to hydrostatic yield point (P*). The combination of hydrostatic overconsolidation and deviatoric overconsolidation tests in this study provides, for the first time, the opportunity to analyze yield curve evolution and compare the effects of purely inelastic volumetric strain (i.e., hydrostatic overconsolidation) against loading with a component of inelastic shear strain (i.e., deviatoric overconsolidation).

Effect of Hydrostatic Vs Deviatoric Strain
As the P* for the deviatoric overconsolidation yield curves (Figure 8b) is not known we cannot normalize the data in the same way as for the hydrostatic overconsolidation tests (Figure 7). However, we can compare curves of similar porosity from the different sets of tests to look at the how the peaks compare. Figure 11a compares the hydrostatically overconsolidated sample at a P* of 70 MPa to the deviatorically overconsolidated sample with an axial strain of 2%. Both these samples have been inelastically compacted to a porosity of 18.0% so should provide a direct comparison between the effects of inelastic hydrostatic versus deviatoric strain.
There is a small difference between the curves with the peak of the deviatorically overconsolidated sample being slightly higher (~5 MPa). However, the difference is minimal and could easily be a result of sample variability. Figure 11b compares the hydrostatically overconsolidated sample at a P* of 130 MPa with the deviatorically overconsolidated sample to an axial strain of 8%. These samples have similar porosities of 14.6% and 14.3%, respectively, so again should provide a useful comparison between the effects of Figure 11. Comparison of similar porosity yield curves deformed by hydrostatic overconsolidation versus deviatoric overconsolidation. (a) The red curve has experienced 2% inelastic deviatoric strain, whereas the black curve has been deformed by hydrostatic compaction. There is a small difference between the curves although this could be as a result of sample variability. (b) The red curve has experienced 8% inelastic deviatoric strain, and the black curve has been deformed by purely hydrostatic compaction. There is a large difference in the yield curve size with the deviatorically overconsolidated curve having a peak that is almost double that of the hydrostatically compacted sample. (c) Schematic diagram to show the evolution of the NCL during deviatoric loading. The hooks (C* 0 ) in the loading curve (red) overlie the initial normal consolidation line (NCL i ) as seen in Figure 5. However, we interpret the normal consolidation line (NCL e ) to evolve during yield curve expansion with deviatoric loading.
Journal of Geophysical Research: Solid Earth 10.1002/2017JB015250 inelastic hydrostatic versus deviatoric strain. In this case, ignoring the change in yield curve shape and appearance of a plateau, the peak of the deviatorically overconsolidated sample is almost double that of the hydrostatically overconsolidated sample. This indicates that there is a drastic effect of inelastic shear strain and anisotropic loading on the evolution of the yield curve, as predicted by Grueschow and Rudnicki (2005).
The evolution observed under deviatoric loading can also help explain the peculiarities observed in the set of experiments from section 2 where the CSL appears to overlie the NCL ( Figure 5). Figure 11b shows the plateau that developed during deviatoric loading sitting directly above the P* of the hydrostatically loaded sample of similar porosity. This is exactly the same phenomenon as that seen in Figure 5 if the CSL is now considered to be a critical-state zone rather than a line. However, in reality, the P* of the deviatorically overconsolidated sample has also evolved as the yield curve had expanded. Therefore, the NCL will migrate as the sample is subject to deviatoric loading, as is schematically shown in Figure 11c. Thus, the peculiarities of Figure 5 are predominantly a result of yield curve evolution in response to the path of inelastic shear strain.
As well as a marked difference in the relative sizes of the yield curves between the hydrostatically and deviatorically consolidated samples, there is also a significant difference in the microstructural evolution. Deformation along the hydrostatic loading trajectory is accommodated by homogeneous compaction (Figure 9b), whereas deviatoric loading leads to a microstructure consisting of multiple sets of localized conjugate shear bands (Figure 9c). Although it is not surprising that the microstructure evolves differently along different loading paths, this must be the reason that such drastic differences in yield curve evolution are observed ( Figure 11). Despite the platy morphology of bassanite making it difficult to determine the exact micromechanical mechanisms in operation during deformation, there is clear evidence in the microstructures (Figure 9) that support differences in mechanical evolution in relation to the loading path. It has been shown previously that different microstructural parameters such as pore geometry and orientation (Bubeck et al., 2017;Griffiths et al., 2017) and the nature of grain contacts (Louis et al., 2009) can produce significant strength anisotropy in porous rock. If microstructural properties evolve significantly during accumulation of inelastic shear strain, as in the deviatoric loading path of our experiments, this could produce a considerable change in the rock strength and in the shape of the yield curve. Models for porous rock deformation typically consider porosity as the only microstructural parameter; however, it is apparent that a much more comprehensive understanding of pore and grain geometry is required to understand fully the nature of yield and subsequent strength evolution.

Potential Implications for Localization
One of the greatest challenges in porous rock deformation is trying to understand the conditions that lead to the onset of localized deformation. In this study we have observed localized shear bands forming in porous bassanite under deviatoric loading conditions (Figures 9c and 12). Many plasticity models for porous rock deformation (e.g., Carroll, 1991;DiMaggio & Sandler, 1971;Schofield & Wroth, 1968) assume an associative flow law, which should only allow localization under strain-softening conditions (Hobbs et al., 1990) that typically accompany dilation (i.e., the low-pressure side of the yield curve). However, this is an oversimplification as localized deformation features have been found to occur under conditions of dilation (Du Bernard et al., 2002), shear (Mair et al., 2000) and also compaction (Mollema & Antonellini, 1996;Olsson & Holcomb, 2000). Localized deformation has also been reported at high pressure, which should suppress instabilities and promote distributed deformation, in studies of both porous (Hirth & Tullis, 1989;Wong et al., 1992; Journal of Geophysical Research: Solid Earth 10.1002/2017JB015250 Zhang, Wong, & Davis, 1990a) and crystalline rocks (Byerlee & Brace, 1969). The inability of associative flow models to predict localization has led to the development of nonassociative flow models, which disassociate the mode of inelastic strain from the yield curve (e.g., Rudnicki & Rice, 1975). Separating the mode of inelastic strain from the yield condition can permit localized deformation to occur under stress conditions where it would not typically be expected in associative flow models. This requires a plastic potential function with additional constitutive parameters that are separate to the yield function in order to describe the nature of inelastic strain at a particular point on the yield curve. Such a framework is outlined for porous rock by Issen and Rudnicki (2000), based on the plastic constitutive parameters of Rudnicki and Rice (1975), to provide a theoretical predication on when the different types of localized feature should occur (i.e., dilation, shear, or compaction band). However, significant discrepancies are found between theoretical predictions and the observed mode of localization in experiment (Baud et al., 2006;Wong et al., 2001). This is attributed to the constitutive model, which is based on nonassociated flow, not being able to capture the different damage mechanisms in operation during porous rock deformation. The complex interplay between pore collapse, grain crushing, and the growth of microcracks will likely be strongly influenced by the local microstructure within the rock such as pore geometry and the nature of grain contacts. The microstructural properties of the rock will in turn affect its localization behavior. For example, it has been previously shown that compaction band formation is highly sensitive to grain size distribution (Cheung et al., 2012). Despite this, it is apparent that we are still lacking a comprehensive understanding of the causes of localization during porous rock deformation.
Although this study does not aim to address specifically whether associative or nonassociative flow is most applicable to the study of porous rock, the recognition of differences in yield curve evolution along different loading paths could have important implications for both, and also for the localization behavior of the material. In the case of associative flow an evolving yield curve might help explain occurrences of localization under conditions that would not typically be predicted from the initial yield criterion. For example, in Figure 11b, the peak of the yield curve for the deviatorically overconsolidated sample overlies the hydrostatic yield point (P*) of a hydrostatically overconsolidated sample of similar porosity. This highlights the variation in yield curve size and shape that can occur along different loading trajectories; which will in turn affect the predicted mode of deformation at a particular point in stress space. The development of a plateau at the peak of the yield curve along the deviatoric loading path could also potentially facilitate localization if compaction and dilation are occurring together as we hypothesized earlier in this section. For example, in the initial experiments in section 2, we observed strain localization and a switch from compaction to dilation in the experiments at the lowest effective pressures of 10, 30, and 50 MPa (Figure 4). This manifested itself as sample-scale shear band with many sets of conjugate shear bands being found in the bulk of the sample away from the main localized feature (Figure 12a). Alternatively in the tests at the highest effective pressures (70 and 90 MPa), the mechanical data suggest that the samples only experienced distributed cataclastic flow. However, even in these samples, less well developed localized shear bands begun to form (Figure 12b), which suggests that dilation was occurring simultaneously with compaction. Perhaps the yield curves of these samples have evolved in a way similar to Figure 10b and the sample now sits in a transitional region of stress space where both compaction and dilation are occurring at the same time. The occurrence of dilation during strain hardening could then help facilitate the formation of localized shear bands. It has also been suggested in a previous study that localized compaction bands also form at the transition region between brittle and ductile deformation in porous rock (Wong et al., 2001).
An evolving yield curve also has important implications for models that assume nonassociated flow. For example, the study of Issen and Rudnicki (2000) uses the internal friction (μ) and dilatancy factor (β) constitutive parameters of Rudnicki and Rice (1975) to predict the type of localized feature that should occur during porous rock deformation. However, the internal friction parameter (μ) is determined from the slope of the yield curve in P-Q space , which is shown in this work to evolve as inelastic strain is accumulated (Figure 8b). In the supporting information, both parameters are calculated as the yield curve evolves along the deviatoric loading path. The internal friction parameter is found to decrease with the formation of a plateau at the peak of the yield curve. Interestingly, the constitutive parameters that were calculated along this loading path always predict that the sample is most likely to localize by forming shear bands, which was observed in all tests where localization occurred (Text S2 and Figure S4). However, irrespective of whether an associated or nonassociated flow law is most applicable, it is apparent that the evolution of Journal of Geophysical Research: Solid Earth 10.1002/2017JB015250 yield curves is important and this is controlled by microstructural evolution in response to different inelastic loading paths.

Conclusions
A series of deformation experiments were performed on porous bassanite to investigate the shape and evolution its yield curve. The curve was found to not be perfectly elliptical in shape as often assumed by critical-state soil mechanics. Rather, the high-pressure compactive side of the curve is partly composed of near vertical limb, which implies that there is a region in P-Q space where the deformation is insensitive to shear stress.
The yield curve shape evolves in response to both inelastic volumetric strain and inelastic shear strain. Samples of the same porosity but with different deformation histories show significant differences in strength. Those that have experienced anisotropic loading with a large component of inelastic shear strain have a yield curve with a peak that is considerably higher than those that have deformed under hydrostatic conditions and accumulated purely volumetric strain. This is a result of a drastic change in the microstructure as the sample begins to localize deformation and form multiple sets of conjugate shear bands in response to inelastic shear strain. A region termed "pressure-insensitive compaction" forms at the peak of the yield curve as the sample deforms along the typical deviatoric stress trajectory of an axisymmetric compression test. In this region, separating the brittle and ductile parts of the yield curve, both compactive and dilational processes occur together, which could help explain previously reported occurrences of localized deformation under conditions that would typically be associated with distributed shear-enhanced compaction. An evolving yield curve is also important to understand for any future models that allow nonassociated flow as this can have implications for the constitutive parameters that describe the inelastic deformation. Future studies should consider the effects of inelastic shear strain on the postyield strength evolution for a range of porous rocks. It is also important to try and constrain what microstructural parameters control the shape and evolution of yield curves and the implications of this for the onset of localized deformation. Our work gives new insight into how previous strain history has a marked influence on porous rock strength, which will in turn influence such processes as the expulsion of aqueous fluids and hydrocarbons in basins, and the movement of fluids in deeper hydrothermal settings.