Model for (De)Compaction and Porosity Waves in Porous Rocks Under Shear Stresses

An understanding of instantaneous and long‐term compaction of porous rocks is important for reservoir engineering and Earth sciences. Experience from laboratory triaxial compression tests and from subsurface operations indicates that shear and volumetric deformations are interdependent. Their mutual dependence results in shear‐enhanced compaction and shear‐induced dilation under short‐term and long‐term loading. Using a classical averaging approach, we consider the evolution of a single fluid‐filled pore in a solid elastoplastic or viscoplastic matrix under combined pressure and shear loading to introduce a new failure envelope and 3‐D constitutive relations for both rate‐dependent and rate‐independent deformation of porous rocks. Our model provides a simple description of rock behavior under a wide range of strain rates. The model predictions agree well with experimental data from triaxial instantaneous and creep tests. Analytical and numerical solutions for solitary porosity wave propagation in viscoplastic rocks in the presence of shear were obtained. New solutions show that new rheological laws have serious implications for porosity waves. Plasticity onset leads to compaction‐decompaction asymmetry and the formation of elongated channel‐like porosity waves. Shear‐induced dilation facilitates porosity wave propagation at fluid pressures below the lithostatic stress. This makes porosity waves a viable mechanism in the formation of focused fluid flow structures in crustal rocks.


Introduction
Many engineering and natural processes in the Earth involve coupled rock deformation and fluid flow (Cai & Bercovici, 2013;Keller et al., 2013;Petrini et al., 2020;Yarushina et al., 2013). Models describing fluid flow in deformable porous rocks can be based in part on the principles of irreversible thermodynamics . Still, the formulation of closure relations represents a significant challenge. The compaction relation that describes volumetric deformation and porosity evolution during external loading is particularly important because it directly affects fluid flow. The combination of irreversible thermodynamics and microscale physics places some constraints on the form of compaction relation . It is well known that the short-term response of most of the compacting porous rocks is elastoplastic (Fortin et al., 2007;Guéguen et al., 2004;Labuz et al., 2018;Makhnenko & Labuz, 2016), while, on a longer time scale, rocks also exhibit viscous properties (Brantut et al., 2013;Chang & Zoback, 2009;Hangx et al., 2010;Makhnenko & Podladchikov, 2018;Renner et al., 2001;Sabitova et al., 2019). Viscous deformation affects the performance of petroleum reservoirs, leading to wellbore stability issues, fault reactivation, and seafloor subsidence (Maranini & Brignoli, 1999). In geotechnical engineering, viscous deformation affects the stability of rock slopes and underground structures, such as tunnels and nuclear waste repositories (Fabre & Pellet, 2006;Ghanbarzadeh et al., 2015;Ma & Daemen, 2006;Tsai et al., 2008). Constitutive models that represent the rate-dependent stress-strain response are implemented in the commercial software (Crook et al., 2008) and are used in numerical studies of reservoir depletion, subsidence, and associated microseismicity (Angus et al., 2015;Yarushina et al., 2017). Viscous deformation and strong coupling between fluid flow and geomechanical deformation can eventually lead to the formation of focused fluid flow, often evidenced in the Earth as dikes, veins, volcanic diatremes, or seismic chimneys (Räss et al., 2014;Minakov et al., 2017). Seismic chimneys are of particular economic importance, as their presence was repeatedly identified as being the most likely indication of leakage and expulsion of gas, oil, and, possibly, CO 2 .
Classical theories of poroelasticity and poroviscosity as formulated by Biot (1962) and Mckenzie (1984) treat shear deformation and compaction as two independent processes. During shear deformation, stress deviator occurrence of critical stress is marked as plasticity onset and the manifold of such yield points is considered a yield surface (Curran & Carroll, 1979;Sevostianov & Kachanov, 2001). Other attempts postulate a specific simplified form of the yield surface and then use ideally plastic solutions to constrain the material parameters used in these laws (Green, 1972). Thus, the evolution of the yield surface during loading, that is, hardening or softening behavior and complete stress-strain relation, cannot be captured by such models. Stress-strain, or compaction, relations are usually fitted using complicated hardening laws involving a number of fitting parameters (Carroll, 1991). The approximate relation between the growth rate of an isolated void and imposed nonhydrostatic stress and strain was developed for simplistic cases of rigid-perfectly plastic materials (Gurson, 1977;Rice & Tracey, 1969) and nonlinear viscous materials (Budiansky et al., 1982). These models are based on the simple problem of the growth (contraction) of an isolated spherical or cylindrical void in an infinite solid matrix. Duva and Hutchinson (1984) use a similar spherical model with approximate solutions to obtain complete stress-strain relations for nonlinearly viscous materials. A viscous damage theory for porous low-cohesion rocks was derived by Ricard and Bercovici (2003), who considered matrix containing isolated voids. This model provides the complete stress-strain relation and predicts enhanced compaction rates in the presence of shear. Triaxial deformation of elastic-viscoplastic porous solid was numerically studied by Koplik and Needleman (1988) on an array of spherical voids. Their results reflect faster compaction rates in the presence of shear.
In this paper, constitutive equations are derived for an incompressible isotropic matrix material containing a dilute distribution of elongated cylindrical pores. The short-and long-term behavior of porous rock is accounted for by considering both rate-independent elastoplastic and rate-dependent viscoplastic deformation. Our RVE consists of a single cylindrical pore subject to spatially homogeneous fluid pressure and nonhydrostatic far-field stresses representing plane strain loading conditions ( Figure 1). We explicitly account for elastoplastic (viscoplastic) deformation of a single pore and its shape evolution during deformation. As this is an already complicated problem, we choose the simplest dilute distribution averaging scheme. This scheme ignores the interaction between the pores and, thus, might be inaccurate at high porosities. However, we show that our model accurately reproduces published experimental data. In the elastoplastic regime, deformation and stresses around the pore can be described using the classical analytical solution of Galin (1946). For the viscoplastic regime, we derived a new analytical solution by modifying elastic and elastoplastic solutions using the viscoelastic correspondence principle. We assume that the initial yield of the matrix material on the microscale is determined by the von Mises criterion. The interested reader can find analytical solutions for the stress and strain fields in the representative volume in Appendix A. Averaging of the analytical solution provides the relation between volumetric strain and stresses for elastic, elastoplastic, and viscoplastic regimes (see Appendix B for complete derivations). The obtained macroscopic compaction relations can be rewritten in a form compatible with the plastic flow theory and, thus, can be further used to derive the functional form of the failure envelope. Because we consider both contraction and dilation of the pores, our solution provides a complete failure envelope that accounts for tensile, shear, and compaction failure.

Macroscale Constitutive Equations for Porous Rocks
Classical viscous models in geodynamics assume that viscous deformation is accumulated in response to any stress over time. The stress-strain rate behavior remains the same for high and low stresses (Gerya, 2019;Mckenzie, 1984;Stevenson & Scott, 1991). On the other hand, existing viscoplastic models that were applied largely to metals assume that viscoplastic flow occurs only when stresses reach yield stress (Chaboche, 2008;Perzyna, 1966;Perzyna & Drabik, 1989). Laboratory experiments show that viscous deformation can be accumulated even at low stresses and that stress-strain rate behavior changes at higher stress levels An infinite body of incompressible elastoplastic (viscoplastic) material containing a long cylindrical pore with radius R is subject to a uniform remote stress σ ∞ x ≠ σ ∞ y , σ ∞ xy ¼0. A uniform pressure p is applied at the pore boundary in fluid-saturated porous material. In dry rocks, the pore boundary is traction-free. (Maranini & Brignoli, 1999;Sabitova et al., 2019;Tsai et al., 2008). Viscous and viscoplastic constitutive equations admit stresses within, on, and outside the yield surface. Plastic constitutive laws admit stresses only within and on the yield surface, which limits the unrealistic stress buildup possible in viscous models. Thus, to keep classical limits and be consistent with experimental observations, we assume that the total strain rates _ ε ij can be split into elastic (i.e., reversible), plastic (i.e., time-independent and irreversible), viscous, and viscoplastic (i.e., time-dependent and irreversible) parts as follows: Thus, the total inelastic deformation is a sum of instantaneous plastic strain and accumulated viscoelastic and viscoplastic strain, reflecting the evolution of the microstructure. Viscoplastic deformation is independent of the loading history and is a function of the current stress state only, while plastic strains depend on the loading history. Volume averaging procedures show that in homogeneous isotropic rocks, elastic volumetric deformation is unaffected by shear stresses (Appendix B, see also Yarushina, Räss, et al., 2015). Thus, the usual incremental poroelastic constitutive equations that are fully consistent with Biot's theory can be used for _ ε e ij : where p ¼ φp f þ 1 − φ ð Þp s is the total pressure, p f ,p s are fluid and solid pressures, φ is the porosity, τ ij is the total stress deviator, _ ε e is the volumetric elastic strain rate, _ γ e ij is the deviatoric elastic strain rate, α BW is the Biot-Willis coefficient, K d is the drained bulk modulus, and G d is the drained shear modulus, while dot denotes the material time derivative (see Table 1 for notations). Positive sign of strain indicates compression. These equations were derived multiple times based on principles of irreversible thermodynamics and effective media theory (e.g., Coussy, 2004;Lopatnikov & Cheng, 2004;Yarushina & Podladchikov, 2015, and references therein). Plastic constitutive equations can be written using the nonassociated flow rule: where Q = Q(p e , τ ij ) is the plastic potential, _ λ is the nonnegative plastic multiplier, and p e ¼p − p f is the effective pressure. The plastic flow rule must be accompanied by the yield criterion and standard plastic consistency conditions, which together limit stress buildup.
Viscous strain rates do not require threshold stress and are linearly related to stresses (e.g., , and references therein), namely, where η b and η s are the bulk and shear viscosities, which, among other things, might depend on temperature and porosity. Our micromechanical model suggests that the simplest porosity dependence has the form where C is the geometrical factor, which is equal to 1 for cylindrical pores, 4/3 for spheres, and ≪1 for elliptical or crack-like pores (Bercovici et al., 2001;Schmeling, 2000;. Often, viscoplastic strain rates are described by means of a viscoplastic flow rule (Chaboche, 2008;Heeres et al., 2002;Perzyna, 1966;Perzyna & Drabik, 1989). Our micromechanical model results in the following simplest viscoplastic flow rule: where Y is the yield stress and W is the viscoplastic potential that will be defined in the next section. The gradient of W determines the direction of viscoplastic flow. A nonnegative viscosity consistency parameter Λ ≥ 0 specifies the magnitude of viscoplastic strain rates. Its values must be calibrated experimentally.

Yield Criterion
In porous rocks, the yield criterion is a somewhat arbitrary notion, as it might correspond to the onset of local plastic failure around stress concentrators (such as in the models of Curran &Carroll, 1979 andSevostianov &Kachanov, 2001), or it might correspond to the point where plasticity is significant and plastic zones around individual imperfections have coalesced (e.g., Green, 1972). In our geometrically simple model, plasticity onset around individual cylindrical pores in a matrix obeying the von Mises or Tresca yield criterion will be defined by the rhomb centered around point p e = 0,τ = 0 (black dotted line in Figure 2a). If smaller imperfections are present, as in dual porosity materials, the matrix will obey the Mohr-Coulomb yield criterion, and plasticity onset around a single pore will be determined by the same rhomb shifted along p e axis by p 0 = C 0 tanϕ, where C 0 is cohesion and ϕ is the angle of internal friction of the matrix material (black dashed triangle in Figure 2a). However, plasticity onset around an individual pore is not representative of the macroscopic yield condition. Yield onset around individual pores first results in a nonlinear deformation. Macroscopic plastic flow occurs when growing plastic zones around individual pores start to significantly influence the stress-strain response of the rock. This is illustrated in triaxial compaction experiments, where rock failure is monitored by acoustic emission. Indeed, acoustic emission starts much earlier than macroscopic failure of the sample exhibits itself by deflection of a compaction curve (Fortin et al., 2006(Fortin et al., , 2009Stanchits et al., 2011). Such deflection is accompanied by a significant increase in acoustic emission. Thus, the shape of the compaction curve provides important indications of the plasticity onset. Let us note that yield function enters into elastoplastic stress-strain relations after the plastic flow rule 3 is complemented with the standard plastic consistency condition (e.g., Chakrabarty, 1987;Yarushina et al., 2010;Gerbault et al., 1998)    pore in a nonhydrostatic stress field (see Appendix B for derivations). The resulting yield function takes the form is the equivalent shear stress. The shape of the corresponding failure envelope is characterized by six independent material parameters: n, m, τ 0 , Y, p 0 , and α τ (Figures 2a and 2b). Linear hardening modulus h defines the evolution of the yield surface, which corresponds to the isotropic expansion of the initial yield surface. Parameters Y and p 0 define tensile and compaction strength. The tensile failure limit is defined by the critical stress Y d = p 0 − Y, while critical pressure for the onset of pore collapse is given by Y c = p 0 +Y. Usually, in experiments, this pressure corresponds to the upsurge in the acoustic emission activity (e.g., Fortin et al., 2009). Because pores and imperfections are present on all scales, decompaction strength in geological materials is different from pore collapse pressure. In the absence of imperfections of smaller size, for example, in porous metals, the failure envelope would be symmetrical with respect to the τ and p axes and α τ = 0, p 0 = 0. Exponents n and m define the curvature of the failure envelope. These depend on the pore geometry and their interaction. In the material with very low porosity composed of aligned cylindrical pores, n = 2, m = 1, and τ 0 ¼Y = ffiffi ffi 5 p . The failure envelope for these values of parameters is shown as a blue line in Figure 2a. Deviations from idealized geometry lead to different values of these parameters and, thus, to different shapes of the failure envelope. Parameter α τ defines the asymmetry of the failure envelope with respect to the p axis (Figure 2b), that is, different plastic moduli in extension and compaction triaxial tests. Note that the derived failure surface encircles rhombs describing failure onset around individual pores.

Viscoplastic Flow Potential
Viscoplastic flow potential may depend on the yield function, W = W(F). In this case, the viscoplastic flow rule is associative and viscoplastic flow is normal to the yield surface. For nonassociative viscoplasticity, W might depend on some other stress function. Here, viscoplastic flow potential W was chosen in a form that both reproduces the theoretical compaction relation derived from the micromechanical model and reduces to zero at the yield surface. The latter requirement guarantees the continuity of strain rates with failure onset. We suggest the following functional form of the viscoplastic potential: Function f is obtained by taking the logarithm of F/Y+1 when hardening is ignored. Thus, it represents an alternative form of the yield function. One can see that f = 0 when F = 0. The shape of surface W = 0 is quite like failure envelopes ( Figure 2c). Substitution of Equations 9 and 10 results in the following associated viscoplastic flow rule. where To avoid singularities that lead to numerical issues, we suggest using even values of m. In the simplest case of m = 2 and α τ = 0 Combining viscous and viscoplastic strain rates, we obtain the following total rate-dependent, or creep, part For the simplest case of m = 2 and α τ = 0 These equations are rewritten in a simpler form that is reminiscent of previous viscous equations used in many analytical and numerical models (Connolly & Podladchikov, 2007;Gerya, 2019;Omlin et al., 2018;Schmeling et al., 2012;. Here, is the effective total bulk viscosity of visco-viscoplastic rocks, is the dilation pressure, and 10.1029/2020JB019683 Journal of Geophysical Research: Solid Earth is the effective shear viscosity. Our model differs from most of the viscoplastic constitutive laws in that it allows for viscous flow at stresses below the critical yield, although with much higher viscosity. This is consistent with experimental data showing that viscous strains would accumulate at both high and low stresses (Crook et al., 2008;Makhnenko & Podladchikov, 2018). Plasticity onset reduces bulk and shear viscosities significantly, leading to enhanced deformation rates (see Figure 3). Shear and bulk strain rates depend on both effective pressure and shear stress. Increasing τ leads to increased compaction rates ( Figure 3a) and further reduction of effective bulk viscosity ( Figure 3b). This is known as shear-enhanced compaction and is observed in creep experiments as well as in the usual triaxial compaction tests (Skurtveit et al., 2013;Xiao et al., 2006;Xiao & Evans, 2003;Zhu et al., 1997). As in other viscoplastic models, more complicated hardening laws and more nonlinearities can be introduced.

Compaction-Dilation Boundary
In our model, the transition from compaction to dilation depends on the interplay between elastic, plastic, and rate-dependent deformation. For rate-dependent deformation, the compaction-dilation boundary can be directly found from Equation 16 by putting _ ε r ¼0. For the simplest considered case of m = 2 and α τ = 0, the compaction-dilation boundary will be defined as As f depends on both p e and τ, this defines a nonlinear curve, which closely follows the failure envelope at p e < Y. It is shown as a thick black curve in Figure 3a. As effective pressure increases further, the compaction-dilation boundary follows a straight line, which will be vertical in the case of α τ = 0 and more inclined at higher values of α τ .
The onset of dilation at positive effective pressures is due to the presence of shear stresses. This is known as shear-induced dilation, which is commonly observed in compaction experiments (e.g., Makhnenko & Labuz, 2015;Vajdova et al., 2004). Several groups have recently performed studies involving multiloading triaxial creep tests on various rock types (e.g., Sabitova et al., 2019;Tsai et al., 2008;Zhang et al., 2016). In these experiments, the axial stress is applied stepwise to cylindrical rock samples while the lateral confinement is held constant. After each load increment, the specimen is allowed to creep for a time interval ranging from hours to days, during which deformations are continuously monitored. The experimental results show that the volumetric strain curves exhibit volumetric compaction at low stresses, followed by dilatancy at higher stresses. The strain rate during compaction is smaller than the strain rate during dilation. Moreover, significant strains were developed during dilation before macroscopic failure of the sample (e.g., Sabitova et al., 2019;Tsai et al., 2008). These results support our model, which predicts viscous dilation at positive effective pressures and leads to compaction/decompaction asymmetry.

Compaction Data
We aim to derive viscoelastoplastic constitutive equations to model the rheology of geological media with applications to fluid flow in mind. However, before we proceed further to fluid flow applications, we compare the derived yield surface and compaction relations with some of the available laboratory results.
In this paper, we do not aim for complete experimental verification of our model, which would require a good dataset with compaction and creep data on the same material under various relevant loading conditions. This lies outside the scope of the present paper. Most of the rock mechanical tests in the literature are performed under triaxial conditions. Very few papers present plane strain results (Makhnenko & Labuz, 2014;Makhnenko & Labuz, 2016). Thus, for comparison, we choose compaction data from triaxial experiments from the literature (Baud et al., 2006;Tsai et al., 2008). In compaction triaxial experiments, a cylindrical sample is initially loaded hydrostatically (σ 1 = σ 2 = σ 3 ) at drained conditions until the desired level of confining pressure p c = σ 2 = σ 3 is reached. After that, p c is kept constant and only axial load σ 1 is further increased to facilitate shear loading so that σ 1 > σ 2 = σ 3 . The mean effective stress and the differential stress measured in experiments correspond to our p e and τ = σ 1 − σ 3 , respectively. Volumetric deformation in compaction experiments is estimated as porosity reduction ε = (φ 0 − φ) Á 100%. In our model, the compaction relation for the triaxial stress conditions takes the form (see Appendix B for derivation)

10.1029/2020JB019683
Journal of Geophysical Research: Solid Earth In triaxial experiments, loading follows a specific trajectory in the plane, so that dτ dp e ¼ 0; Substituting Equations 8 and 27 into Equation 26, one obtains the following compaction relation, which might be used directly for the description of experimental results: In the hydrostatic limit, this it reduces to These compaction equations and this failure envelope were implemented into MATLAB code, which is used further for comparison of theoretical results with experimental data.
The experimental data that we use for comparison with model output were taken from Baud et al. (2006). They describe the compaction of Darley Dale sandstone with an initial porosity of 13% (see Figure 2a in (Baud et al., 2000)). Hydrostatic compaction was performed at effective pressures up to 500 MPa. Five nonhydrostatic curves were obtained by initial hydrostatic loading up to 150, 180, 200, 240, and 300 MPa, with confining pressure, p c , kept constant at these levels ( Figure 4). At low effective mean stress, data show a nonlinear rapid compaction, which might be explained by the closure of microcracks. A crack porosity can be introduced here to capture this portion of a compaction curve (Shapiro, 2003). For simplicity, we do not model an exact compaction trend below the crack closure pressure. Instead, we approximate this stage with almost instant porosity reduction, assuming that the elastic compaction of stiff porosity starts after ε = 0.75%. After crack closure pressure, hydrostatic compaction data show almost linear dependence of volumetric Journal of Geophysical Research: Solid Earth deformation on effective pressure until pore collapse slightly enhances the compaction rate, leading to nonlinearity beyond critical stress p * = 380MPa (black stars in Figure 4). The presence of shear stresses reduces the critical stress for the onset of pore collapse and significantly influences the slope of the compaction curves. Compaction data from Baud et al. (2006) are supplemented by experimental failure envelopes shown in Figure 4b with an asterisk (ε p = 0), circles (ε p = 0.1%), and diamonds (ε p = 0.25%) for different levels of accumulated volumetric plastic strain. The best fit for all experimental curves was obtained using G = 1.4 GPa, Y = 190 MPa, p 0 = 190 MPa, τ 0 = 128 MPa, n = 5, m = 1, α τ = 0.2, and h = 12.5 GPa. Figure 4 shows theoretical predictions (solid lines) plotted on top of the respective data points for six compaction curves corresponding to different levels of confining pressure and three failure envelopes corresponding to three different levels of accumulated plastic strain. Note that once fracturing or strain localization is initiated, full 2-D or 3-D numerical simulations of elastoplastic deformation and elastic unloading are required to describe the experimental results (Stefanov et al., 2011).

Viscous Creep Data
Viscous, or time-dependent, deformation of porous rocks is studied in a number of triaxial creep experiments recently summarized in a review paper from Brantut et al. (2013). These studies show that time-dependent deformation can be very significant in shales, rock salt, tuff, sandstones, and limestones but also present in hard rocks such as granite and basalt. Most of the creep experimental studies focused on stress-strain-time behavior and the influence of temperature and a pore fluid chemistry on creep. Some studies also report axial creep strain rates as a function of applied stresses (Brantut et al., 2013;Heap et al., 2015;Rybacki et al., 2015). The works of Zhang et al. (2015), Makhnenko and Podladchikov (2018), and Sabitova et al. (2019) report volumetric strain rates. Relatively few studies focus on characterizing the yield surface and viscoplastic potentials (Maranini & Brignoli, 1999;Tsai et al., 2008;Weng et al., 2010).
In most of the materials, development of the creep strain in time is characterized by three distinct stages: the primary, or transient, creep phase with strain rate decreasing with time; the secondary, or stationary, creep phase with a constant (i.e., time-independent) strain rate; and a final, accelerating stage followed by specimen failure. Tsai et al. (2008) performed multistage triaxial creep experiments on dry samples of Mushan sandstone from Taiwan, in which both volumetric strain and shear strain were measured. Experiments were performed under a constant room temperature and the time of each experiment was adjusted so that the stationary creep stage was reached. Elastic deformation was carefully subtracted from the total strain by using loading-unloading-reloading cycles. Vectors of plastic and viscoplastic strain increments were obtained during the short-term and creep stages of experiments, respectively, and plotted on the (p e , τ)-plane. After that, plastic and viscoplastic potential surfaces were plotted, assuming that they are normal to the direction of the strain increments. Tsai et al. (2008) showed that obtained viscoplastic potential has the same functional form as the plastic potential and yield surface. Maranini and Brignoli (1999) also documented the viscoplastic failure envelope during creep tests on porous limestone, which is very similar to the standard failure envelope.
Here, we compare our viscoplastic flow potentials with those obtained by (Tsai et al., 2008). Figure 5 shows a comparison of the experimental data with the theoretical viscoplastic potential defined by Equation 9 with n = 2, m = 2, and α τ = 0.2. For the green curve, we used p 0 = 37MPa, Y = 21MPa, and τ 0 = 33MPa; for the blue curve, we used p 0 = 43MPa, Y = 27MPa, and τ 0 = 46MPa; and for the red curve, we used p 0 = 46MPa, Y = 35MPa, and τ 0 = 57MPa. The results show reasonable agreement.

Porosity Waves in the Presence of Shear
Porosity waves were suggested to be a fast and efficient fluid transport mechanism forming various focused fluid flow structures (Appold & Nunn, 2002;Barcilon & Richter, 1986;Cai & Bercovici, 2013;Jordan et al., 2018;Mckenzie, 1984;Räss et al., 2014;Richard et al., 2012;Scott & Stevenson, 1984;. Most of the previous models ignored the influence of deviatoric stresses on wave propagation, assuming that there is no shear or deviatoric stresses in the rock. Yet, recent numerical models that include full mechanical formulation with all components of the stress tensor show that propagation of porosity waves is also associated with changes in shear stresses, especially ahead of the wave head (Omlin et al., 2018). However, these models used ad hoc rheology with decompaction weakening only. Here, we investigate how derived rheology influences the propagation of porosity waves. First, we consider the 1-D steady-state compacting porous column and derive an analytical solution for purely viscous and viscoplastic porosity waves in the presence of shear. After that, we look at the transient 2-D numerical solution for porosity waves in a medium, which viscoplastic behavior is described by using the constitutive equations derived above.

1-D Steady-State Viscous Solution With Shear
We consider equilibrium compaction or uniform sedimentation of a vertical column of fluid-saturated porous rock. During sedimentation, low-density fluid segregates from the deformable porous matrix. Over time, compaction in the upper part of the sedimentary column evolves toward a steady state, while the porosity profile becomes independent of time and forms a stationary porosity wave (Figure 6a). One-dimensional stationary solutions can be obtained by putting ∂/∂t = 0 and replacing material time derivatives with d/dt = v∂/ ∂z, where v = v(z) is a vertical component of solid velocity. Deviatoric stresses are present in the media. In our derivation of the one-dimensional steady-state solution, we follow the same procedures as in previous publications from Connolly and Podladchikov (2000) and Yarushina, Podladchikov, and Connolly (2015), which provide them in more detail. Here, we offer a brief summary of the governing equations and show how they are modified in the presence of shear.
Mass conservation for fluid and solid phases gives the following nontrivial components of the solid velocity and Darcy flux: where φ b is the constant background porosity, q s is the sedimentation rate, and q 0 is the background Darcy where τ z is the vertical deviatoric stress, g is the gravity, and ρ ¼ ρ f φ þ ρ s 1 − φ ð Þis the total density. The only nontrivial component of the shear stresses can be related to the vertical velocity using the linear viscous law: and, thus, the total force balance becomes At the remote boundary, we assume the lithostatic pressure gradient, that is, p ′ z ¼ ∞ ð Þ¼ − gρ. Equation 33 must be complemented by the Darcy law where k is permeability, μ is fluid shear viscosity, and p f and ρ f are the pressure and density of the fluid. Introducing effective pressure p e and using Equation 30, the Darcy law may be written in the form We assume that, at the remote boundary, dp e /dz = 0; thus, from Equation 35, it follows that where k b is the constant background permeability and Δ ρ = ρ s − ρ f is the difference between solid and fluid densities. Combining Equations 35 and 33, we arrive at the following hydraulic equation: In the limit of d 2 v/dz 2 → 0, it coincides with the hydraulic equation from the models of Connolly and Podladchikov (2007) and Yarushina, Podladchikov, and Connolly (2015). This equation must be solved together with Equation 30 for solid velocity and the equation for viscous compaction: Assuming power law porosity-permeability dependence Equations 37-39 can be combined into a single ordinary differential equation, the so-called phase portrait equation:

10.1029/2020JB019683
Journal of Geophysical Research: Solid Earth is the characteristic pressure scale, is a nondimensional speed of the porosity wave, and is the ratio of shear and bulk viscosities. Its solution is given by the following expression: Equation 44 coincides with the solution previously given by Yarushina, Podladchikov, and Connolly (2015) if α = 0. Depth profiles of porosity and pressure may be obtained numerically by solving coupled ordinary differential equations where

10.1029/2020JB019683
Journal of Geophysical Research: Solid Earth is the compaction length. Figure 7 shows the porosity wave solution to Equations 44-46 for α = 100, φ b = 0.01, υ = 3, and ω = 4. The comparison of viscous solutions with shear (α > 0) and viscous solutions without shear (α = 0) in Figures 7 and 8 shows that shear stresses reduce pressure generated within the porosity wave. The nonlinear dependence of permeability on porosity plays an important role in this reduction, especially for rocks with a strong porosity dependence of permeability (i.e., with high values of υ). For example, shales generally have large permeability exponents, in the range υ = 24.9 − 55.5 (Dong et al., 2010). David et al. (1994) report υ = 4.6 − 25.4 for five different sandstones. Our results show that for υ = 3, pressure reduces about two times in the presence of shear (Figure 8a) in comparison to a reference case without shear (Figure 8c). For υ = 20, pressure is reduced over 30 times in the presence of shear (Figure 8b). Similarly, the amplitude of the wave, that is, the maximum porosity that the wave can reach, is affected by υ and the magnitude of the shear stresses. Porosity inside the wave is much higher for rocks with low permeability exponent υ, while in rocks with higher υ, porosity inside the wave does not increase as much. For example, wave amplitude φ/φ b ≈ 13 for υ = 3, while it is only φ/φ b ≈ 1.3 for υ = 20. Yet, implications for permeability change inside the wave can be quite significant, as even a small porosity change inside the wave can increase permeability orders of magnitude at high values of υ in Equation 39. It is interesting to note that the presence of shear stresses may also lead to sharp pressure gradients inside the wave (Figure 8a) when the difference between permeability exponent υ and nondimensional wave speed ω is too high. These pressure gradients at the pressure front can be responsible for numerical instabilities in the numerical codes based on poroviscoelastic models with full stress tensor, as implemented by Omlin et al. (2018) and Räss et al. (2019). Instabilities in these codes arise at high values of the permeability exponent, which brings strong nonlinearity into the system of equations. Our 1-D calculations for various values of υ and ω presented in Figure 8 indicate that it is not high values of υ in themselves that cause such problems but, rather, the difference between υ and ω. For example, the pressure profile for υ = 20 looks smoother (Figure 8b) than the pressure profile for υ = 3 (Figure 8a) for ω = 24.

1-D Steady-State Viscoplastic Solution With Shear
Now we will investigate how plastic deformation alters porosity wave propagation. As in the previous section, we will look at the upward rising porosity wave in a setup shown in Figure 6a. The only driving force is buoyancy arising due to the fluid being less dense than the rock. All previous porosity wave solutions assumed explicitly or implicitly that p e = 0 at the remote boundary (e.g., Connolly & Podladchikov, 2007;Rice, 1992). We relax this assumption by allowing fluid pressures to be smaller than lithostatic and assuming, instead, that p e = p d > 0 at the remote boundary, that is, effective pressure corresponds to the dilation pressure in our viscoplastic constitutive law given by Equations 20 and 21. This assumption is still somewhat idealistic; in the future, other boundary conditions must be tested.
Rising fluid pressure at the top of the wave will induce viscoplastic failure leading to dilation, even at positive effective pressures. At the lower part of the wave, pressures are reduced and, thus, will fall into the purely viscous domain. This will lead to compaction-decompaction asymmetry, where different constitutive laws will describe the compacting and dilating parts of the wave. Most of the derivations of the previous section will be valid for the viscoplastic case. In the compaction domain, that is, at p e > p d , the viscous porosity Equation 38 from the previous section will still hold, and pressure and porosity evolution will be given by Equations 40, 45, and 46 from the previous section. With a new nontrivial pressure boundary condition p e (φ b ) = p d at the remote boundary, the solution to Equation 40 will take the form In the dilating part of the wave, the viscoplastic compaction Equation 20 might be used to obtain the porosity equation of the form 10.1029/2020JB019683

Journal of Geophysical Research: Solid Earth
For simplicity, we assume that η eff and p d are constant. The first assumption means that in Equation 5, changes in porosity are ignored. The second assumption means that the compaction/dilation boundary is approximated with a vertical straight line in the plastic domain (see Figure 3a), that is, that an elaborate expression for p d in Equation 23 is approximated with a constant value. This is reasonable for α τ = 0 in Equation 8. In the next section, we present a numerical solution that includes the full dependency on pressure and shear stress. Here, we want to explore how the introduction of plasticity and dilatancy at positive pressures affects the porosity waves. Equation 37 does not include assumptions about the rheology of the rock and can be used to derive the hydraulic equation. Substituting Equation 49 into Equations 30 and 37 gives The porosity Equation 49 can be rewritten in the following nondimensional form: Excluding z from these equations, we obtain Integration of the last equation provides an analytical solution for the phase diagrams at p e ≤ p d , namely, The viscoplastic porosity wave solution is presented in Figure 9. It shows that the onset of dilatancy at critical pressure p d significantly alters the porosity and pressure in comparison with a purely viscous solution. As decompaction plastic viscosity is much smaller than purely viscous compaction viscosity (Figure 3), the wave becomes asymmetric, which in 2-D leads to the generation of elongated channels, as we show in the next section. Most importantly, the effective pressures associated with the porosity wave propagation remain positive, that is, generation of a porosity wave in viscoplastic rocks is possible at fluid pressures below the lithostatic stress.

2-D Viscoplastic Numerical Solution
We implemented derived rate-dependent constitutive equations into 2-D numerical code for porous flow in deformable porous viscoplastic rocks using the model formulation from . Our numerical solution ignores elastic deformation. We used the same numerical scheme as described in (Omlin et al., 2018;Poliakov et al., 1993;Räss et al., 2019). The code is available online (at https://doi.pangaea.de/10.1594/PANGAEA.909658). As an initial setup, we use a rectangular domain of viscoplastic deformable porous rock, containing an elliptical high-porosity inclusion in the lower part of the domain formed by, for example, partial melting in magmatic rocks, petroleum generation in source rocks, or dehydration reactions in metamorphic rocks (Figure 6b). Inside the inclusion, the initial porosity is three times higher than the background value, φ b = 0.05. Note that unless inertial (e.g., wave propagation) or elastic processes are involved, fluid flow problems in viscoplastic media do not require initial conditions for pressures. Thus, we do not prescribe any specific values for either fluid or effective pressures. The size of the 10.1029/2020JB019683

Journal of Geophysical Research: Solid Earth
computational domain is 20 × 40 with the unit of compaction length. We impose free-slip boundary conditions for the solid velocities (no shear stress and zero normal velocity). For the fluid flow problem, we assign no flux conditions on the lateral boundaries and constant inflow (bottom) and outflow (top). The only external force acting on the domain is gravity pointing downward. The parameters used in our 2-D simulations are summarized in Table 2.
To compare the results of our model with those of previous studies, we conduct numerical simulations with both viscoplastic and viscous bilinear rheology. Figure 10a shows channels formed in bilinear viscous rocks with a compaction viscosity 100 times higher than the decompaction viscosity. The 2-D simulation results show that our new viscoplastic rheology also leads to the formation of elongated tubular channels (Figure 10b), which are very similar to channels formed in bilinear viscous rocks with decompaction weakening, where the ratio of compaction/decompaction viscosities is R = 100. The rheological asymmetry in the viscoplastic case is caused by a viscoplastic failure at elevated fluid pressures (low effective pressures). It is interesting to note that no negative effective pressure is observed during channel propagation in viscoplastic rocks (Figure 10b). This is related to the new viscoplastic constitutive relations given by Equation 20, which cause dilation of the rock at p e = p d , reducing pressure buildup above dilation pressure p d . The results of 2-D numerical simulations also confirm predictions of the 1-D steady-state analytical solution presented in Figure 9, which shows that the effective pressures associated with porosity wave remains positive. This is different from purely viscous

10.1029/2020JB019683
Journal of Geophysical Research: Solid Earth rocks, in which channel propagation requires fluid pressures higher than the total pressure (compare Figures 10a and 10b). Thus, the new rheology leads to the formation of channelized fluid flow without introducing negative effective pressure, as was the case in previous models. Comparing both models at the same time step, we see that the fluid channels are slightly bigger and rise faster in the viscoplastic model than in the viscous bilinear model. This implies that the new viscoplastic rheology promotes fluid channel formation more efficiently than a decompaction weakening rheology from previous works with R = 100. Another interesting observation is that shear stresses are unevenly distributed during wave propagation in both viscous bilinear and viscoplastic case and, thus, the assumption for the hydrostatic stress state can be very inaccurate in numerical codes.

Journal of Geophysical Research: Solid Earth
Five models with different grid resolutions have been run for the sensitivity analysis. The evolution of maximum porosity in the channel (Figure 11a) and the position of the wave front (Figure 11b) show quite similar trends in all the models. The maximum porosity at the initial time steps demonstrates the major difference among the models. This is due to the strong fluid flow within the reservoir at the beginning of the simulation. The lower the resolution is, the bigger the error it introduces. However, little difference can be seen from grid resolution 192 × 384 to 256 × 512. Figures 11c and 11d illustrate the effect of resolution on the porosity evolution. At nondimensional time = 0.2, two major fluid channels form in all three models, and the wave front tends to be higher in the low-resolution model than in the high-resolution model. The runs with low resolution diverge from the correct trend. Three small fluid channels can be seen in the low-resolution model, while only two small channels are found in the other two models. Similar observations were reported for porosity wave solutions by Räss et al. (2019). Thus, considering both the statistical values and model images, our model achieved good convergence when the grid resolution increased.

Discussion and Conclusions
An essential feature of our constitutive model is that it builds on the effective media theory. We used an analytical solution for the deformation of a single pore in a nonhydrostatic stress field in our derivation of fully 3-D viscoelastoplastic constitutive laws. Yield surface, flow potential, and compaction relations are obtained by simple volume averaging of an analytical solution. Parameters in the compaction relations are generalized to account for more complex pore geometry. Our model accounts for the (1) stress dependency of material parameters; (2) weakening of the rock with the onset of failure; and (3)  Porosity waves were previously proposed as a mechanism for generating fast and focused fluid flow in the shallow and deep Earth, which is often evidenced in the form of dikes, veins, pockmarks, mud volcanoes, fluid escape pipes, or gas-conducting chimneys, among other structures. The contrast in the compaction and decompaction viscosities is an important parameter that controls the mode of fluid flow focusing on viscously deforming porous rocks (Connolly & Podladchikov, 2007). In rocks with the same compaction and decompaction viscosity, fluid flow instability gives rise to porosity waves in the form of spherical blobs. Such blobs can propagate upward as a self-sustained body at a speed that exceeds the background porous fluid flow rate but that is still too inefficient to explain the formation of focused fluid flow structures (Barcilon & Richter, 1986;Mckenzie, 1984;Schmeling, 2000). In rocks with different compaction and decompaction viscosities, fluid flow instability leads to the formation of elongated chimney-like structures (Connolly & Podladchikov, 2007;Omlin et al., 2018;Räss et al., 2014;. The ratio between compaction and decompaction viscosity determines the length of the chimney and the speed of its upward propagation. However, compaction/decompaction asymmetry in the previous models was achieved in an ad hoc manner by postulating constant but different viscosities at positive and negative effective pressures. Additionally, these models predicted that fluid pressure must exceed lithostatic pressure to produce a porosity wave. These deficiencies limited the success of these models. Our new model also has different compaction and decompaction viscosities. However, the difference is achieved naturally as a result of developing plastic failure, as we showed in the previous sections (see Figure 3). Because the transition from compaction to dilation in our model happens at positive effective pressures, fluid pressure does not need to exceed lithostatic pressure to generate a porosity wave. Thus, in our model, porosity waves are a viable mechanism for the generation of focused fluid flow structures. Porosity waves can even compete with hydraulic fracturing, which is often assumed to be responsible for the formation of focused fluid flow (e.g., Arntsen et al., 2007;Bodvarsson et al., 2003). Another application in which the proposed model might be useful is the modeling of shear localization and nonlinear changes in seismic velocity associated with failure processes. This we leave for future work.

Appendix A: Nonhydrostatic Elasto(visco)plastic Deformation of Cylindrical Pores
Often, the choice of an RVE is dictated by available analytical solutions. For compaction and decompaction of porous rock under combined pressure and shear loads, cylindrical pore embedded in an infinite incompressible elastoplastic (or viscoplastic) solid that is subject to homogeneous pressure on the pore wall and a couple of far-field compressive or extensional forces ( σ ∞ x ≠ σ ∞ y , σ ∞ xy ¼0 ) can be considered as an RVE (Figure 1). In the elastic regime, deformation of the pore can be obtained using Kolossof-Muskhelishvili's method of complex potentials, as shown in (Yarushina & Podladchikov, 2007). In the elastoplastic regime, pore collapse is described by the classical analytical solution of Galin (1946), which is extensively discussed in the literature (Chakrabarty, 1987;Kachanov, 1971;Yarushina et al., 2010). Both solutions can be modified to account for the viscous deformation of the solid using the viscoelastic correspondence principle. Changes in pore size and geometry are particularly important for constraining (de)compaction. Thus, here we account only for results reflecting displacements or velocities around the pore. An interested reader can find the full solution with stress distribution in the references above.

A.1 Elastic (Viscous) Deformation of the Pore
In the elastic regime, displacements around the cylindrical pore take the form (Yarushina & Podladchikov, 2007) where G is the elastic shear modulus and are the total pressure prescribed at the far-field boundary, shear stress, and effective pressure, respectively. Due to the mathematical similarity of elastic and viscous problems, the viscous solution can be obtained from Equations A1 and A2 by replacing displacements u r , u θ with velocities v r , v θ and elastic shear modulus G with shear viscosity η (Goodier, 1936).

A.2 Elastoplastic (Viscoplastic) Deformation of the Pore
With increasing load, the plastic region develops around the pore (Figure 1). This corresponds to grain crushing and pore collapse in laboratory experiments. This moment marks the onset of shear-enhanced compaction. We assume that material around the pore obeys Tresca or von Mises yield criteria, both of which take the form where Y is the yield stress in pure shear. Plastic yielding begins at the pore boundary when external loads reach critical condition 2 τ j j þ p e j j ¼ Y : When the plastic region fully engulfs the cavity, stresses and displacements in an elastoplastic rock are given by Galin's analytical solution. Here, we use the similarity of pure elastic and pure viscous formulations to modify the classical elastoplastic solution and find stresses and velocities around the cavity in a viscoplastic rock. We consider separately elastic (viscous) and plastic regions.

A.2.1 Solution for the Plastic Zone
Following Galin (1946), we assume that plastic flow in the vicinity of the pore has radial symmetry. In this case, three unknown stress components must satisfy two force balance equations and the yield criterion, which means that stresses can be found without involving any rheological laws. Therefore, stress distribution in the plastic zone is the same for elastoplastic and viscoplastic materials and is given by the following modification of the original Galin's equations: where the new plastic limit Journal of Geophysical Research: Solid Earth ξ = − 1 for decompaction. However, no analytical solution in a closed form for displacements in the plastic zone is available for σ ∞ x ≠ σ ∞ y .

A.2.2 Solution for the Elastic (Viscous) Zone
The solution in the elastic zone is found using Kolossof-Muskhelishvili's method combined with the use of conformal mapping so that (Yarushina et al., 2010) 2G where ς is a complex variable and Due to the mathematical similarity of viscous and elastic problems, Equations A12-A16 can also be used in the viscous zone if displacements and elastic modulus are replaced with viscous counterparts.
The elastoplastic (viscoplastic) solution is valid as long as shear stresses are relatively small (Yarushina et al., 2010). At high effective pressures p e , the imposed shear load should not exceed approximately 40% of the yield strength (τ=Y < ffiffi ffi 2 p − 1). At low effective pressures (approximately), the admissible shear load grows exponentially with scaled effective pressure that guarantees that the cavity is completely enclosed in the plastic boundary:

A.2.3 Estimation of the Kinematic Field Within the Plastic Region
Though velocities and displacements in the plastic zone are unknown, some estimations of the relative pore contraction (expansion) can be made if the solid rock is incompressible. The conservation of mass in the solid volume of an RVE requires that where Σ is the boundary of the solid volume V s with an outward pointing normal n, ρ is the constant density, and v is the velocity vector. For the two-dimensional problem considered here, Equation A18 can be rewritten as follows: where S s is the area occupied by the solid in the (x,y) plane. When the solid contains a cavity, the line integral in A19 is a sum of the integral over the pore boundary L v and the integral over the external boundary L ext The first of these integrals represents the flux through the pore with area S v and also gives changes in the area of the pore dS v /dt. The second integral is the flux through the total area S. Equation A20 essentially states that in the incompressible RVE with a single pore, all volume (area) changes are due to pore contraction (expansion). Furthermore, in the incompressible solid, any volume containing the cavity would change by the same amount. In particular, where S p is the area enclosed in the elastoplastic (viscoplastic) boundary L (Figure 1). Displacements (velocities for viscoplastic case) at the plastic boundary L can be found from Equations A12 after substitution of complex potentials A13-A15: To express u x , u y through coordinates in the physical z plane, we replaced terms containing 1/ς and 1/ς 2 using relations . The solution to Equation A15 with respect to ς gives where plus and minus signs describe two different half-planes. On the plastic boundary and, therefore, the square root in Equation A24 can be replaced using After this operation, Equation A22 on the plastic boundary may be rewritten in the form To evaluate integral A19, we notice that Substitution of A27 into A28 after some algebra leads to dθ: Taking a contour integral of expression A29 over one quarter of the ellipse gives 10.1029/2020JB019683 Journal of Geophysical Research: Solid Earth ∮ Lp − u y dx þ u x dy ¼ 4 ∫ Lp=4 − u y dx þ u x dy ¼ − πc 2 Y 2 þ 5τ 2 GY ξ : According to A21, dS v /dt can be obtained by taking time derivatives of both parts in Equation A30 for elastic rocks, leading to Appendix B: Averaging of the Solution For a chosen model of RVE, the degree of compaction can be estimated via changes in the pore volume fraction as follows :

B.1 Elastic Compaction
For the circular cavity For elastic rocks, differentiation of radial displacement A1 with respect to time and subsequent substitution into B2 gives where dot stands for the time derivative. Taking into account that S v = πR 2 , the porosity equation for poroelastic rocks becomes G φ dφ dt ¼ − dp e dt : Equation B4 implies that, for small elastic deformations around cylindrical pores, shear deformation does not affect porosity. Effective compressibility in poroelastic rocks is a function of porosity only:

B.2 Viscous Compaction
A viscous compaction relation can be obtained by the substitution of u r with v r and G with shear viscosity η s in Equation B3, resulting in As in the elastic case, the presence of shear stresses during viscous deformation does not influence compaction. Effective viscosity takes the form

10.1029/2020JB019683
Journal of Geophysical Research: Solid Earth

B.3 Nonhydrostatic Elastoplastic Compaction
For elastoplastic deformation, dS v /dt is estimated using Equation A31, which, together with S v = πR 2 , gives the elastoplastic porosity equation in the form Upon substitution of Equation A16 for c, the elastoplastic compaction relation may be rewritten as with effective elastoplastic compressibility Equation B9 reduces to elastic porosity Equation B4 when there is no shear stress and effective pressure reaches the critical value |p e | = Y for plasticity onset.

B.4 Nonhydrostatic Viscoplastic Compaction
Replacing displacements and elastic moduli with their viscous counterparts in Equation A30 for the rate of the pore volume change and substituting it into Equation B1, one obtains Substitution of c with Equation A16 gives the final form of the viscoplastic porosity equation in the presence of combined pressure and shear loading: This equation determines effective viscosity of the form In the absence of shear, it reproduces the effective viscosity given in our previous work .

B.5 Derivation of Yield Criterion
The compaction Equation B9 may be used to derive yield function. For this, we will rewrite elastoplastic constitutive equations for plane strain conditions: Plastic multiplier _ λ can be found from the consistency condition ∂F ∂p e dp e þ ∂F ∂τ dτ þ ∂F ∂ε p dε p ¼ 0 Direct comparison of Equations B17 and B9, and subsequent integration of Equation B17 with respect to F, gives yield function of the form which gives derived compaction equations if τ 0 ¼Y = ffiffi ffi 5 p , n = 2, m = 1, α τ = 0, p 0 = 0, and h = G. Indeed, substituting Equation B19 into Equation B17, we obtain