Thermodynamic Characterization of Chemical Damage in Variably Saturated Water-Active Shales

A constitutive framework is developed for variably saturated water-active swelling rocks undergoing chemical damage using modified mixture theory and continuum damage mechanics. The Helmholtzian thermodynamic potential for the skeletal system is derived as a function of the state variables including deformation, damage, two-phase fluid pressures, and chemical potential. Using this, in addition to chemo-poroelastic constitutive equations, a thermodynamically consistent first-order estimation of the damage variable is developed. The working of the theory is shown through the numerical example of water uptake in clay-rich shale rocks solved by the finite element method. The numerical results portray the significance of including variably saturated conditions in constitutive equations as a unique damage-dependent poroelastic behavior was observed for wet and dry regions. The theoretical-based damage estimation corroborated by previous experimental observations illustrates that the rock strength is dominantly controlled by the time of exposure to water rather than the level of water saturation. Contrary to what was perceived, the results show that poroelastic and chemo-poroelastic responses do not coincide even in less reactive shales due to the time-dependent water-induced microstructural deterioration of the rock. The microstructural deterioration increases the storage and flow capacity in the water-saturated region giving rise to substantive spatio-temporal changes in matrix stresses. The research findings provide valuable insights to understand how poromechanics plays a role in causing water uptake in water-sensitive rocks and how such behavior is coupled with associated microstructural chemical damage. •Helmholtzian thermodynamic potential is used to derive damage evolution in shales. •Rock strength controlled by time of water exposure rather than level of water saturation. •Poroelastic and chemo-poroelastic responses differ significantly even in less reactive shales. •Damage increases flow and storage causing significant stress response changes. •Damage occurs only in water-saturated region, and rock tortuosity delays damage.


Effective stress in non-damaged space L ijkl
Fourth-order elastic stiffness tensor kl Second-order deformation tensor F Helmholtz free energy W Helmholtzian thermodynamic potential p w/g Water/gas pore pressure Chemical potential of specie Y D Damage conjugate variable v w Water pore volume v g Gas pore volume m bound Bound mass of specie

Introduction
Subsurface porous structures can be saturated with one or more fluids depending on the depth of burial and pressure-temperature conditions (Siddiqui et al. 2014;Ahmed 2018). Coupling between the flow of such fluids and mechanical deformation of the porous structure is critical in subsurface events such as hydrocarbon recovery, hydrology, surface subsidence, seismicity, nuclear waste disposal, carbon dioxide sequestration, etc. (Jha and Juanes 2014). One of the most widely used coupling schemes between fluid flow and geomechanics is the Biot theory of poroelasticity (Biot 1941). This theory initially dealt only with porous materials saturated with a single fluid with later extensions to account for chemical, thermal, and plasticity effects for single fluid-saturated porous materials (Heidug and Wong 1996;Ghassemi and Diek 2003;Roshan and Fahad 2012;Roshan and Oeser 2012;White et al. 2017). However, the challenge persisting is related to simple two-phase poroelasticity and its extension to other effects evident from a dearth of literature (Borja 2006;Chen 2013a;Song and Borja 2014;Choo et al. 2016). Generalization of the theory is important considering the usual existence of two-phase mixtures (gas/water, water/oil, gas/oil) in deeper rocks for geothermal energy extraction, CO 2 geological sequestration, and hydrocarbon production. Also, for shallow rocks close to the atmosphere, pore spaces are often partially saturated, i.e., filled by air and water where a general two-phase poroelasticity theory can be useful. From a modeling standpoint, constitutive equations for poroelasticity are often used in numerical simulations of coupled flow and geomechanics problems. To develop such constitutive theories, the concept of mixture theories is used. Mixture theory fundamentally describes the interaction of the constituents of a mixture (Truesdell and Toupin 1960;Green and Naghdi 1965;Goodman and Cowin 1972;Passman 1977). In the jargon of continuum mechanics, a mixture exists when each point is co-occupied by two or more constituents (Katsube and Carroll 1987). Based on this assumption the conventional mixture theories exhibit challenges when applied to fluid-saturated porous media because when fluid is removed, a solid material is left instead of a dry porous solid skeleton (Katsube and Carroll 1987). A modification was, thus, introduced to solve this problem. In, what is now referred to as, the modified mixture theory, dry and fluid-saturated porous solid skeletons were replaced by an equivalent homogenous solid and an equivalent homogenous liquid (Katsube and Carroll 1987). The equivalent solid and liquid act as two unique continua occupying the same point. The approach is now commonly used (Coussy 2005) but is not always referred to as the modified mixture theory. Heidug and Wong (1996) were among some early researchers who derived, using the modified mixture theory, the single-phase chemo-poroelastic theory based on non-equilibrium thermodynamics and continuum mechanics for porous geomaterials. Their method was picked up and followed by many more constitutive modeling studies. In contrast to the classical mixture theory (Truesdell and Toupin 1960;Truesdell and Bowen 1984), the method of Heidug and Wong (1996) does not distinguish explicitly between solid and fluid phases. This approach was different from the consolidation theories of Terzaghi (1943) and Biot (1941) where many complex constitutive models for realworld applications including unsaturated porous media were based on (Meroi et al. 1995;Sanavia et al. 2002;Li et al. 2006;Seetharam et al. 2007). The issue with consolidation theories is that extending the constitutive equations to incorporate complex chemical fluid-solid interactions becomes challenging (Laloui et al. 2003;Chen and Hicks 2013). This drawback is easily overcome by the modified mixture theory approach (Heidug and Wong 1996). For example, Chen and Hicks (2010) extended the concept of Heidug and Wong (1996) to develop constitutive equations for an unsaturated rock with chemical aspects added to the modified mixture theory (Chen and Hicks 2013).
Chemical effects represent an internal coupling phenomenon where the force transfer is not apparent. Although microlevel interactions could be captured with molecular dynamics modeling, such an approach is not feasible to extend the coupling effects to macro-scale continuum models (Chen and Hicks 2013). Here, by macro-scale, it is implied the scale at which the hydro-mechanical coupling occurs. The modified mixture theory (Heidug and Wong 1996) uses continuum mechanics along with non-equilibrium thermodynamics to establish a link between mechanical and chemical coupling by looking into energy dissipation mechanisms and the time evolution of state variables (Heidug and Wong 1996;Chen and Hicks 2010). This approach was used lately to extend the unsaturated poroelastic constitutive equations to incorporate chemical (hydrational swelling, osmosis) and thermal effects (Chen 2013b;Chen et al. , 2016Chen et al. , 2018a. The drawback of these studies and others (Zienkiewicz et al. 1990;Meroi et al. 1995) is that for unsaturated porous media, the gas phase (air) pressure was ignored, i.e., average pore pressure was used for an averaged fluid mixture and chemical-induced damage (microstructural deterioration of the rock fabric) was not considered which can be significant in clay-rich rocks such as shales. Thus, such an approach again may not apply to real-world engineering problems. The issue with assuming gas phase pressure to be atmospheric is that it implies the gas phase can flow freely without any resistance. In fluid flow simulations, the results obtained without such simplification have been found to differ considerably from the assumption of a passive gas phase (Morel-Seytoux and Billica 1985). This can be significant where a gas phase (especially for gases other than air) is present which can have a considerable change in flow and water swelling properties of medium with changes in pressure and/or temperature. Moreover, neglection of chemical damage implies ignoring the permeability variation that occurs due to the evolution of micro-cracks brought about by osmosis or clay swelling. Therefore, a rigorous modeling framework where gas flow and chemical damage are considered is important especially for clay-rich rocks. Schrefler and Scotta (2001) made progress with regards to considering an active gas phase where two continuity equations for the water and gas phases were used. Efficient numerical stability was portrayed in their study using the various fluid relationships between saturation, capillary pressure, and relative permeability. More recently, Cheng (2020) developed phenomenological constitutive equations for unsaturated porous media and presented simplified observable relationships between fluid saturation (in the case of partial saturation) and laboratory-measurable constants (in the case of full saturation). For pure mechanical effects, consistency was found between the phenomenological micromechanical approach and the thermodynamics-based variational energy approach. The solution to the constitutive equations was not presented but the developed model seemed robust. Although such an approach is appealing towards efficient numerical and engineering applications, it is not easy to extend it to couple other multi-physics phenomena including thermal, or reaction mechanisms such as chemical damage (Cheng 2020).
Chemical damage is a critical multi-physics phenomenon that induces the microstructural deterioration of the water-active clay-rich rocks e.g., swelling induced damage in shales due to clay-water reactivity. Rocks, like any solid material, are composed of a unique microstructure. The mechanical behavior of the material such as strength, stiffness, etc. depends on this microstructure (Wei et al. 1989).
Observations have shown that a change in microstructure results in the reduction of the mechanical properties (e.g., strength) and can eventually lead to cracking (failure) of the material (Panteleev et al. 2021). When this microstructural change (e.g., micro-cracks or micro-voids) is due to the clay-water reactivity in shales, we refer to this as the chemical damage in this study. The generated micro-cracks need not necessarily always cause measurable permanent deformation of the material but at the same time can lead to a significant loss in strength of the material even without the presence of any mechanical load (i.e., elasto-damage or brittle damage (Khan et al. 2004(Khan et al. , 2007Murakami 2012;Browning et al. 2017;Jia et al. 2020)). It is the accumulation, growth, and localization of this damage that causes the final material to fail. At the scale of the representative elementary volume (REV), a macroscopic continuous damage field is used to represent the microstructural defects. We refer the readers to the pioneering works of continuum damage mechanics for further details (Kachanov 1958;Lemaitre 1985Lemaitre , 1992Lyakhovsky et al. 1997b;Murakami 2012).
From a thermodynamics point of view, the principle of local state asserts that a point in a continuum can be described by specific state variables. As mentioned, damage refers to the deterioration of the material's microstructure. This microstructural deterioration at a point in the continuum is an internal state and can be represented by an internal variable called the damage variable, denoted by D (Lemaitre 1992; Murakami 2012) which is a scalar assuming the microstructural discontinuities are isotropic. What inspires us to thermodynamically include chemical damage in the chemo-poroelastic constitutive theory is the specific focus of this study-water loss in shale matrix. Shale rock is a porous medium characterized by pores and a solid skeleton where the pores are initially filled with gas (with irreducible water). After hydraulic fracturing, there are main hydraulic fractures and several hydraulically induced micro-fractures where the water resides. During the well shut-in period, water imbibes from the main fractures and the induced micro-fractures into the gas-saturated matrix. Since the induced micro-fractures are hydraulically poorly connected to the main fracture (Dahi-Taleghani and Olson 2011;Nandlal and Weijermars 2019;Siddiqui et al. 2021), the imbibition modes from the main fracture and induced micro-fractures are different. From the induced micro-fractures, spontaneous imbibition is the main mechanism of water uptake by the shale matrix (Tokunaga 2020). The extent of this spontaneous imbibition in causing water loss in shales is the focus of the investigation in this study. As the water imbibes, two-phase flow occurs in the shale matrix and micro-cracks develop in the rock due to water-shale (clay minerals) interactions. These micro-cracks are a timely result of the chemically driven microstructural alteration of the matrix and should not be confused with the above-mentioned hydraulically induced micro-fractures generated during hydraulic fracturing.
It is clear that though advancements have been made to the thermodynamics-based rederivation of the Biot poroelasticity theory to complex porous media scenarios including two-phase flow conditions (e.g., Coussy 2005), attention is still warranted to develop a thermodynamically consistent theory explaining the chemical damage and its association to water loss problem in water-active shales (Roshan et al. 2015;Siddiqui et al. 2018Siddiqui et al. , 2019. The development of such a theory is the focus of this study. Therefore, in this study, using non-equilibrium thermodynamics and continuum mechanics principles, based on the modified mixture theory, we develop a two-phase chemo-poroelastic, damage constitutive theory. The theory builds on previous constitutive theories to provide a new first-order estimate of shale matrix damage due to water uptake in variably saturated conditions. This can help explain the extent of water loss in the shale matrix and serve as a basis for the development of further advanced constitutive theories. The solution to the developed theory was numerically performed to specifically address the shale water loss problem in hydraulic fracturing operations. The details of the derivations are presented in the next section.

Theoretical Development Using Thermodynamics
When the first and second laws of thermodynamics are applied to any system, they involve the state variables that characterize the system's internal energy. Hence, when the evolution of such a system is of interest, an appropriate framework is provided by these laws to develop the respective constitutive equations. In this study, we apply the thermodynamic laws to a porous continuum to account for different possible couplings in the development of the constitutive equations. Such a thermodynamics-based approach takes inspiration from the pioneering works of Biot (Biot 1955(Biot , 1977Biot and Temple 1972) and Sherwood (Sherwood 1993, 1994, and the comprehensive works of Coussy (Coussy 2005(Coussy , 2010 and Heidug and Wong (1996). The extension to the basic poroelasticity theory is performed in this study by accounting for the second fluid phase, clay swelling, and chemical-induced microstructural deterioration, in the overall Helmholtz free energy (F) of the porous solid. In other words, a Helmholtzian-based thermodynamically consistent constitutive theory for a two-phase chemoporoelastic (hydro-chemo-mechanical) gas-water system is developed. For this, we consider an isotropic macroscopic region of volume V at isothermal conditions within the rock matrix having a boundary R which is allowed to deform only infinitesimally. This is henceforth referred to as the representative elementary volume (REV). The REV contains both solid (grains) and voids. Such a continuum-scale representation of the rock microstructure inherently implies that V is also significantly larger than the largest grain size of the REV. The REV boundary (R) is considered closed to the flow of solid matter while fluid (and/or solute) inflow and outflow can occur. This further implies that R is attached to the external solid (grains) part of V such that the deformation of the solid matrix is governed by the evolution of R in space (Heidug and Wong 1996). It is also considered that the water phase consists of only one type of solute where water and solute exhibit chemical potential w and c , respectively (the associated electrolytic behavior is insignificant). It is also imposed that the gas phase is not adsorbed to the matrix but exists as a free phase.
With regards to the microstructural deterioration of the REV, we reiterate that the focus of the theory is reserved for chemical damage. The readers are advised that mechanical damage may become critical during cyclic stress perturbations which is not of concern in the shale water loss context. Furthermore, the water-clay interaction and consequent swelling are considered to be instantaneous to the time frame of the movement of fluids in the rock. This maintains the local physio-chemical equilibrium and the mechanical equilibrium (Heidug and Wong 1996). It is important to note that the theory is developed for strictly monotonic loading (fluid withdrawal from underground reservoirs), i.e., there is no cyclic loading considered nor are there any drying-wetting cycles. To incorporate cyclic stress (or effective stress) perturbations, thermodynamical consistency dictates new internal state variables related to plasticity and viscous effects will need to be added to the free energy potential of the rock skeleton (Arson 2020;Jacquey and Regenauer-Lieb 2021). Moreover, mechanical damage accumulation with cyclic stress perturbation becomes significant and an additional mechanical damage variable needs to be also added (Kuhl et al. 2000). This will give rise to complex damage coupling between chemical and mechanical damage variables that need extensive experimental characterization (Le Bellégo et al. 2003;Comi et al. 2014).
The state of any system can be represented by a fundamental function called thermodynamic potential. Starting from the first and second laws of thermodynamics, the expression for the differential form of the Helmholtz thermodynamic potential for the rock (REV) skeleton can be derived (see Appendix 1 for the detailed derivation). As discussed in the Introduction, the REV-scale (continuum-scale) constitutive modeling is based on the modified mixture theory wherein the porous medium is the superimposition of the skeleton continuum and the fluid continuum. However, the thermodynamics-based approach usually focuses on the skeleton deformation as it is the one that can be observed (e.g., Heidug and Wong 1996;Coussy 2005). Hence, similar to Heidug and Wong (1996), and Coussy (2005), we too focus on the skeleton system starting with its Helmholtzian thermodynamic potential: The variables , D, p w , p g , constitute the state variables characterizing the state of the skeleton system and they, respectively, represent strain, damage variable, water pore pressure, gas pore pressure, and chemical potential of different species (water, gas, solute). According to the postulate of local state, these state variables are macroscopic variables without referring to the microscopic level (Pokorska 2008). We consider the fundamentals of the energy-based approach to constitutive modeling to be well documented. Therefore, even though a major part of the derivation leading to Eq. (1) is presented in Appendix 1, for a more rigorous step-by-step guide, the readers are referred to other sources (c.f. Appendix 1 in Heidug and Wong (1996), or Chapter 3 in Coussy (2005)).
Taking the time derivative of W in terms of the state variables, we associate the state variables with their conjugate thermodynamic state variables , Y D , v w , v g , m bound : where, from Eq. (2), we can write the state equations as (Algazlan et al. 2022): where is stress, Y D is damage conjugate variable, v w/g is variation in pore water/gas content, and m bound is the mass concentration of species in the skeleton system (e.g., in (1) W = W , D, p w , p g , . (2) any occluded porosity). The bold font represents tensor or vector quantities as per the norms of continuum mechanics. Henceforth, indicial notations (without bold font) will be used, e.g., ij is the second-order stress tensor and ij is the corresponding second-order strain tensor; both based on a Cartesian coordinate system (i, j = 1, 2, 3) . Hence, Eq. (2) can be written as: At this point, the readers are advised that some materials undergo healing (of the accumulated damage) through different mechanical or chemical processes. A thermodynamically consistent theory of damage and healing mechanics was first proposed by Barbero et al. (Barbero et al. 2005;Arson 2020). For thermodynamically consistent healing mechanics, healing is considered a dissipative process that is characterized by an extra internal variable in addition to the damage variable D . We chose not to incorporate healing in our model as we develop the theory for fluids withdrawal from underground reservoirs where stress reversals to aid healing are rare. Healing mechanics is a work in progress and the readers are referred to the comprehensive works of the group of Arson et al. (Zhu and Arson 2016;Shen and Arson 2018; Arson 2020) for details.
Next, differentiating Eq. (3) with respect to the time, the following constitutive equations for the evolution (time rates) of stress, variation content of water phase, variation content of gas phase, variation in interlayer (bound) water, and damage evolution are obtained: where the first terms on the right side of Eqs. (6)-(9) are double dot products. The set of constitutive equations (Eqs. (5)- (9)) also holds at equilibrium conditions. The local state postulate imposes that during any evolution of a system, it is characterized by the same state variables that characterize the equilibrium states (Coussy 2005). In other words, this implies that Eq.
(2) holds even without time dependency. The set of constitutive equations (Eqs. (5)-(9)) are already geometrically linear due to the infinitesimal strain and Cauchy stress assumptions. These equations can be further physically linearized with the understanding that the thermodynamic response coefficients, , and Z are material-dependent constants (Heidug and Wong 1996). For isotropic materials, M ij , N ij , and S ij are diagonal tensors and can be written in terms of scalars w , g , and , respectively: Here, w and g are the Biot coefficients for water and gas, respectively, is the swelling coefficient of each component, and ij is the Kronecker delta. For binary solutions = 0 M c ∕RT where 0 is the common swelling coefficient representing both solute and solvent, M c is the molar mass of solute, R is the universal gas constant, and T is the temperature (Ghassemi and Diek 2003). It is further emphasized that ∑ i=w,g i = , where is the Biot coefficient of the saturated porous medium (Jha and Juanes 2014). It is also noted that i are proportional to the respective saturation of each individual phase ( i = S i ) where ∑ i=w,g S i = 1 (Jha and Juanes 2014) such that in the case of complete saturation with respect to any fluid we will have i = enabling a smooth transition from partially saturated to fully saturated conditions.
The thermodynamic response coefficient, L ijkl then possesses the meaning of the elastic stiffness matrix which is the fourth-order tensor written, for a damaged material, using the strain equivalence theory proposed by Lemaitre (1992): where G is the rock's shear modulus, K its bulk modulus, and D is the damage variable. According to the equivalent strain concept proposed by Lemaitre (1992), the virgin material's elastic constitutive equation can be used for the damaged material by multiplying the elasticity matrix with the (8) (11) damage variable. It is noted that the stiffness tensor formulation is chosen to be in terms of the shear modulus ( G ) as it was observed that under isotropic stresses no microstructural alteration (i.e., damage) occurs (Siddiqui et al. 2020a). Using Eqs. (10)-(11), Eqs. (5)- (7) can be written as: When the saturation of the wetting phase (water) increases to attain full saturation, physically there is no second phase present p g = 0 . In such a case, the constitutive theory transitions to a fully saturated flow of a single wetting phase and the constitutive equation for gas (Eq. 14) becomes redundant and all terms associated with the variable p g vanish in other constitutive equations. The attention of the readers is sought here in realizing that Eqs. (12)-(13) simplify to the classical Biot poroelasticity theory when the presence of gas phase ̇g and chemical reaction ̇ are neglected. Moreover, neglecting only the chemical term ̇ in Eqs. (12)- (14), results in the same unsaturated poroelastic constitutive equations derived by Coussy (c.f. Section 7.3 in (Coussy 2010)). The approach to deriving thermodynamics-based constitutive equations in this study follows Heidug and Wong (1996), which is slightly different from the approach of Coussy (2005Coussy ( , 2010. Nonetheless, the fact that similar consistent constitutive equations emerge from both approaches provides confidence in the accuracy of the developed theory. It is noted that since gas does not cause any swelling, Π will remain zero. Furthermore, as was investigated in the authors' previous study (Siddiqui et al. 2020a), the effect of water/gas pore pressures is negligible in damage evolution, i.e., W = X = 0 and invoking the strain equivalence concept V ij = 0 . Hence, we have:

Field Equations
Before deriving the field equations, the chemical potential in the constitutive equations (Eqs. (12)- (14)) is replaced by postulating an ideal solution (dilute solutions) with a single solute where the solute chemical potential can be estimated as Diek 2002, 2003): where R is the gas constant, T is the system temperature, M c is the solute molar mass, and C c is the mass fraction of the solute which is related to the water mass fraction (C w = 1 − C c ) . The linear approximation of Eq. (16) yields Diek 2002, 2003): where C c is the average solute mass fraction over the range of interest. Rearranging the Gibbs-Duhem equation, the chemical potential of water is then expanded (Heidug and Wong 1996;Chen 2013a;Chen et al. 2018a): where w and c are the density of water and solute relative to the unit fluid volume, respectively.

Momentum Balance for Solid Phase
To arrive at the final field equation for the solid phase, we consider the mechanical equilibrium condition in the absence of the body, gravity, and dynamic forces (Coussy 2005): and incorporating Eqs. (17) where υ is the rock's Poisson's ratio and the solid displacements u i are related to the strain through: (

Transport Laws and Mass Balance for Water, Gas, and Solute
It is recalled that during the theoretical development (Appendix 1), it was stipulated that no interfacial exchange or mixing occurs between water-gas phases. The dissipation function (Eq. (94)) is then used to derive the phenomenological equations with the addition of the thermodynamic driving forces and corresponding fluxes using Onsager's theory of irreversible processes (Onsager 1931). The 'flow' of damage is not included in the formulation of transport phenomenological equations. Such phenomenological equations relate the thermodynamic driving forces to their consequent fluxes (Heidug and Wong 1996;Ghassemi and Diek 2003)-in this case, the water, gas, and solute fluxes. Hence, the phenomenological equations including gas transport can be written as: Here, the thermodynamic driving forces are ∇p w , ∇p g , and ∇( c − w ) . It is commonly considered that the crosspermeability coefficients can be neglected L 12 = L 21 = 0 implying that the pressure gradient of a fluid phase can only generate a flux of that phase. The readers might be interested to read the view of others (Kalaydjian 1990;Mannseth 1991;Avraam and Payatakes 1999) who argue that the cross-permeability coefficients are nonnegligible even though their magnitudes are much lower than the diagonal permeability coefficients L 11 , L 22 . The significance of the cross-permeability coefficients has been widely debated (Avraam and Payatakes 1995). The cross-permeability coefficients represent the viscous coupling exerted between fluid phases which is a strong function of the size of the fluid/fluid interface and the ratio of viscosities of the two fluids. With the pore sizes commonly found in shales (Villamor Lora et al. 2016) and a wide contrast in water/gas viscosities, we ignore L 12 , L 21 in this study. It is also obvious that there will be generally no solute present in the gas phase; hence L 23 = 0 . The effect of water and gas pressure on solute transport is also minimal allowing L 31 = L 32 = 0 (Roshan and Rahman 2013; Marbach and Bocquet 2019). It is noted that Onsager's symmetry . (Onsager 1931) is not inherently invoked here as is mostly done in the literature. Whether the symmetry holds or not is subject to specifically designed experiments which will not be discussed here for now. It is also noted that for dilute solutions which can be considered to be ideal, the following approximation is valid (Ghassemi and Diek 2003): where C c and C w are the average mass fraction of solute and water, respectively. Therefore, the transport equations for water, gas, and solute phases are derived as: where u w , u g , and J c are the water, gas, and solute flux vectors, respectively, D e is the effective solute diffusion coefficient, k is the absolute permeability of the porous media, k rw and k rg are the water and gas relative permeability, respectively. Also, ℜ is the reflection coefficient or the membrane efficiency of the shale which varies from 0 to 1 and w and g are the dynamic viscosities of water and gas, respectively. The relative permeability to be used in the transport laws (Eqs. (27)-(28)) is estimated from the fractal analysis technique which has proven to be efficient in clay-rich rocks (Zhang et al. 2018b, a;Siddiqui et al. 2020b): where S w is the water saturation. The readers are advised that there are other familiar relative permeability models such as van Genuchten (1980), Brooks andCorey (1966), Purcell (1949), etc. which may be used if applicable to the specific shale under consideration. A good comparative view of these and other models can be found in Li and Horne (2006).
The time evolution of variation of compressible/incompressible fluid contents ( ) is expressed as (Ghassemi and Diek 2003): where i = w, g and K f ,i=w,g is the fluid bulk modulus for water and gas, respectively. For incompressible water, as is obvious, the second term on the right of Eq. (33) becomes negligible due to its large fluid bulk modulus. Now, invoking the continuity equation (mass balance for water, gas, and solute) yields: where u w , u g are the fluxes of pore water and gas, respectively. Assuming water is the wetting phase, the capillary pressure, p c is defined as: where Substituting Eqs. (13) and (14) into Eq. (33), and using B = 1 K ( − 1) proposed by Heidug and Wong (1996) ( K here is the bulk modulus of porous medium), the time evolution of variations in water and gas contents are obtained as: Substituting the water and gas fluxes (Eqs. (27) and (28)), and the time change of water and gas contents (Eqs. (39)-(40)), into the water and gas continuity equation (Eqs. (34)-(35)), yields the coupled water and gas two-phase field equations: (39) Finally, the coupling coefficients, M ∕ may have various forms but for now can be defined as (Cheng 2020): where K c = − dS w dp c is the inverse capillary pressure derivative defined using the fractal concepts (Siddiqui et al. 2020b): The solute transfer equation can be obtained by substituting Eq. (29) into Eq. (36). Keeping the convective term c u w ∇C c in Eq. (36), we arrive at:

Damage Evolution
Next, to close the set of field equations, we need to characterize the complementary damage force evolution law (Eq. 15). This requires the knowledge of the thermodynamic response coefficient (Y) . Hence, we alternatively use the entropy production (dissipation) due to microstructural deterioration to directly characterize damage evolution Ḋ . For chemical damage due to the combined action of water-solute-clay interactions, we can formulate the evolution of D , based on the stress constitutive equation (Eq. 12), by defining stress ( ) in − space instead of − space: Similar to Eq. (105) (see Appendix 2 for details), in − space, we have the following expression for the dissipated energy due to an infinitesimal increase of chemical damage at constant strain and fluid pressures: Integrating Eq. (48) and invoking the definition of c ( c = 0 M c ∕RT ), and the Gibbs-Duhem relationship (Eq. 17), the following expression for Y D is obtained after substitution in Eq. 107 (see Appendix 2 for details on Y D which is equated to −Y D to work with positive quantities): Next, the entropy production (dissipation) due to microstructural deterioration is equated to the dissipation by chemical interactions (Amiri and Modarres 2014) to arrive at: Solving for Ḋ using the Gibbs-Duhem relationship (Eq. 17) on the right-hand side for c gives the following expression: where RT converts 0 to molar density units. Hence,Eqs. 20,41,42,46, and 51 form a set of constitutive equations that describe the dynamic evolution of thermodynamic state variables of bulk displacement, gas, and water pressure, solute transfer, and damage in a porous system saturated with water and gas. The set of material constants and other inputs needed to solve the field equations were obtained through the authors' previous experimental studies or the literature (more details in Sect. 3). The readers are advised here that the direct entropy-based chemical damage characterization relies on the time-dependent increase in continuum damage to weaken the material (Amiri and Modarres 2014). The actual force causing damage ( Y in Eq. 15) is not measured which requires in situ monitoring and advanced sensors (Amiri and Modarres 2014). Such an approach is considered practically sufficient for metals, though much research is warranted for other materials.
The damage (micro-cracks) initiation due to clay-fluid interactions is reflected in the absolute permeability on the continuum scale. Based on the experimental observations from the authors' previous study (Siddiqui et al. 2020a); chemically induced micro-fractures only appear considerably in presence of anisotropic stresses. It is therefore assumed that the permeability alteration by these microfractures occurs only during deviatoric loading or shearing. This alteration is given by the commonly used Eq. (52) below (Choinska et al. 2007;Zhang et al. 2007;Zhu and Wei 2011;Qu et al. 2019), which is considered only when anisotropic stresses exist: where D is a constant, k 0 and 0 are the initial permeability and porosity, respectively. Hence, with the evolution of D , the absolute permeability of the matrix (REV) is updated according to Eq. (52). The readers are advised that nonlinear auxiliary relationships like Eq. 52 and the like (e.g., Eqs. 30, 31, 44, 45) do not necessarily cause internal consistency problems in the context of engineering applications such as in this study and others (Laloui et al. 2003;Gens et al. 2006;Cheng 2020). It is fair to state here that in the absence of chemical damage (D = 0) , there will be no change in the permeability of the system. It is noted that the basis for Eq. (52) are experimental observations of permeability and porosity increase or decrease with damage accumulation or healing, respectively (Tenthorey et al. 2003;Lyakhovsky and Hamiel 2007). Such an expression is physically motivated by the logarithmic response of the static friction coefficient with time (Lyakhovsky and Hamiel 2007). The damage evolution measured with acoustic emissions has also been shown to exponentially increase with time (Lyakhovsky et al. 1997a). The porosity ( ) after microstructural alteration is given by (Peng et al. 2015;Qu et al. 2019): where the Δ symbol represents the difference between the current and initial value of a parameter, v is the total poroelastic volumetric strain, and swell is the swelling strain. It is assumed that the swelling strain is recoverable. The readers are advised that different forms of Eq. 53 can be found in the literature (Peng et al. 2015(Peng et al. , 2018, but in this study, the matrix compressive strain component is included in Δ v (c.f. (Aghighi et al. 2021)). The swelling strain to be used in Eq. (53) is derived from the stress constitutive equation by dividing the swelling stress component by the bulk modulus ( K): The formulation of swelling strain, like the chemical damage variable, is a function of the coupled chemical reaction of water and solute, i.e., f (S w , C c ) . In Eq. (54), it is intuitive to say that the contribution of water saturation to swelling is dominant, i.e., even if the imbibing water has lower salinity compared to in situ water, due to the increasing S w with imbibition, there will still be swelling. Although both water and solute actions are responsible for clay swelling, it would be ideal to have their unique damage contributions through having separation functions, i.e., f 1 (S w ) + f 2 (C c ) . This, however, requires further experimental characterizations that are a focus of our future studies. The following section pursues the numerical solution of the developed constitutive equations using the finite element method based COMSOL Multiphysics platform.

Model Description and Numerical Example
Based on the authors' previous studies (Siddiqui et al. 2020a), where it was observed from micro-scale investigations that the deformation induced by fluid exposure was (53) localized and considerably trivial compared to the sample volume, one simplifying observation is made, i.e., the crosscoupling terms ẇkk and ġkk are negligible in Eqs. (41) and (42), respectively. The final field equations implemented in COMSOL are Eqs. 20, 41, 42, and 46, and the damage evolution law Eq. 51, which form the complete set of twophase chemo-poroelastic, damage constitutive theory. These coupled equations are solved through the utilization of different modules and physics interfaces. Details about the different modules and interfaces utilized in this study can be found in COMSOL Users Guide (COMSOL Multiphysics 2009).

Model Conceptualization for Shale Water Loss
The numerical example is presented as a specific case of water loss in shale rock matrix during hydraulic fracturing operations. Representative mechanical loads are applied to the REV. Considering the normal faulting stress regime, as the common case for shale reservoirs, the mechanical loads are the overburden stress and minimum horizontal in situ stress. These are assumed to be the highest and lowest stresses acting on the REV, respectively. Although the intermediate stress is relatively important on the mechanical response of the rocks, the minimum and maximum stress effects are shown to be dominating (Cao et al. 2016;Yu et al. 2020). It is considered that hydraulically induced micro-fractures have long lengths compared to their widths leading to the valid assumption of the plane strain problem. With these arguments, a 2D REV model can simulate the conditions of the shale matrix during well shut-in. The 2D domain geometry (REV) is, thus, shown in Fig. 1.

Initial and Boundary Conditions
The initial water phase pressure in the REV is zero (0) MPa as it exists initially at irreducible water saturation (0.17). The initial gas phase pressure is set as 34.5 MPa (5000 psi) based on the reservoir pressure with a saturation of 0.83. The gas reservoir pressure can vary with the burial depth and can sometimes be exceptionally high due to stratigraphic conditions. The initial and later water saturations are calculated using the auxiliary equations presented earlier where the capillary pressure is known from the authors' previous study [estimated using fractal theory (Siddiqui et al. 2020b)]. The initial gas saturation is obtained by subtracting water saturation from unity. The initial solute concentration of the REV is set as 0.5 mol/m 3 which is dependent on the composition of in situ pore water. The flow boundary conditions for the 2D domain are set as a zero-flux top, bottom, and right boundaries, and constant pressure left boundary (Fig. 1). The constant pressure on the left boundary is the water inlet pressure and the gas outlet pressure. The water pressure in the induced microfractures after well shut-in is usually below the reservoir pressure (Rutqvist et al. 2018) representing spontaneous imbibition conditions. In this study, it is set as 20.7 MPa (3000 psi). Spontaneous imbibition is a slow process and progresses during the period of shut-in, e.g., more than a week (Rutqvist et al. 2018). The imbibed water causes the gas to flow out of the induced micro-fracture-matrix face, i.e., a counter-current imbibition is modeled (Al-Ameri et al. 2018). The solute boundary condition on the left boundary is set as 0.1 mol/m 3 such that the salinity of the in situ water is higher than that of the fracturing fluid. The mechanical boundary conditions are set as a roller for the bottom and right boundaries whereas the left and top boundaries can deform (Fig. 1). It is noted that the value of in situ stresses depends on the burial depth of the gas reservoir and the stratigraphic location. As per sign convention, in the simulation results presented, compressive stresses are negative.

Finite Element Discretization
A finite element mesh is created using mapped quadrilateral elements. The entire 2D geometry is divided into 5600 quad elements with an element-area ratio of 0.1. The mesh was refined near the shale-fracture interface to accurately capture the imbibed waterfronts at early times (Fig. 1). The scatter/line plots of numerical solutions presented subsequently are plotted along the yellow horizontal line aligned with the X-direction of the REV as shown in Fig. 1. The simulation time was 500 h which is around 20 days, within the typical well shut-in durations post hydraulic fracturing. The numerical implementation in COMSOL is summarized in Fig. 2. It is noted that the nonlinear expressions (e.g., Eqs. 30,31,44,45,52) were defined as local user-defined variables in the COMSOL Multiphysics platform. The Constant (Newton) nonlinear solver method of COMSOL is proven to be robust to handle the nonlinearities (Gudala and Govindarajan 2021;Wang et al. 2021;Xue et al. 2021).

Simulation of the Full Constitutive Theory
The validation of the developed constitutive theory is presented in Appendix 3 through a comparison of simplified sub-models with analytical or previously published numerical results. Here, the full constitutive theory is solved and presented for the case of water loss in shale matrix REV. The required inputs are tabulated below (Table 1). The numerical results portray the complete simulation of the developed constitutive theory for two-phase, damage chemo-poroelasticity. The important sensitivity analysis is discussed for parameters that affect damage and also how damage affects the stress response in the rock.
Before discussing the focal issue of chemical damage, we first present the results for water-gas flow with associated solute transport as a part of the two-phase chemo-poroelastic model developed in this study. As expected, the pressure diffusion and consequent waterfront movement into the matrix are slow due to the low permeability of the shale (Fig. 3a,  c). After 500 h (of shut-in), the water imbibition length is around 2 mm or 13.3% of the REV length (Fig. 3c). In the low-velocity spontaneous imbibition process, no dispersion or mixing occurs when water imbibes into a gas-saturated porous medium (Gonfiantini et al. 1998) and hence the water flow follows piston-like displacement.
The chemical effects are incorporated using the solute transport equation. It is these chemical effects that initiate the chemical damage. It is seen that the time to equalize the solute concentrations in the matrix and fracturing fluid (external water) is very long (Fig. 3b). Even after 500 h, the solute outflow reached only < 2 mm. It is also seen that outward solute transport follows closely with water saturation (compare Fig. 3b, c). This is aligned with the developed constitutive theory where only the imbibed water favors solute transport and there is no solute transport in the gas-saturated region of the matrix.
The swelling stress swell generated in the matrix is shown in Fig. 3d. As seen, the maximum swelling stress is around 0.42 MPa (59.5 psi) corresponding to the maximum amount of water imbibed. With the progress of water imbibition, the region of the matrix experiencing chemical swelling expands (compare Fig. 3c, d). It is worth noting that the generated swelling stress is too low to overcome the Fig. 1 Schematic of the 2D geometry, mesh, and boundary conditions used for simulating spontaneous imbibition from induced micro-fractures into shale matrix high magnitudes of confining stresses. Hence, a macro-scale impact of such stresses is not manifested through, for example, the creation of large micro-fractures with significant plastic strains. This was indeed seen in micro-scale experiments in the authors' earlier study (Siddiqui et al. 2020a). No swelling stress was generated in the gas-saturated region since the water is still at irreducible water saturation.
With water and solute transport occurring in the shale matrix, chemical damage is prominent which affects the stress state inside the matrix. The significance of the developed theory is highlighted through the observation of a unique damage poroelastic response which is challenging to observe with other two-phase poroelastic formulations (Chen 2013b;Chen et al. 2018b). The stress distribution  depends on the evolution of damage in addition to the level and extent of water saturation giving rise to unique poroelastic responses for water and gas-saturated regions (dashed lines in Fig. 4). The unique poroelastic response in waterand gas-saturated regions of the rock is evident. There is a clear transition between the two regions which is characterized by what is referred to in this study as 'poroelastic response interface' (PRI). The PRI diffuses further into the matrix with time. As the water-gas interface moves further into the matrix, the PRI also expands. When damage is neglected (solid lines in Fig. 4), the poroelastic behavior is significantly different. It is seen that while the stress in the gas-saturated region remains unchanged (compare dashed and solid lines in Fig. 4), the stresses in the water-saturated region change significantly due to neglecting the damage-induced permeability evolution in the water-saturated region. When damage is neglected, the effective stresses in the damaged region are grossly misestimated (compare the difference between dashed and solid lines in Fig. 4). The difference in stress due to damage and without damage increases with time as the water-saturated region experiences increasing microstructural deterioration, i.e., the rock is weakened due to chemical damage evolution. For example, the difference  David et al. (2015) where mechanical instability (including damage visualized by X-ray computed tomography) was found to be associated only with regions that imbibed water. This experimental observation corroborates the developed theory.
The damage variable profiles give insight into the evolution of damage in the water-saturated region. The damage variable evolution is independent of the confining conditions. It is seen that after 5 h of water imbibition, the chemical damage variable D = D chem reached a maximum of 0.13 at the fracture-matrix interface and diffusively reduces further into the matrix before neutralizing in the gas-saturated region (Fig. 5a). The value of D chem is always highest at the fracture-matrix interface due to the maximum water saturation levels therein. At 50 h, D chem increased to 0.28 and by 500 h has already reached the theoretical limit of 0.3 amounting to localization of microstructural alteration. The permeability in the water-saturated increases significantly with time due to damage evolution while the permeability in the gas-saturated region remains constant (Fig. 5b). Permeability is seen to increase ~ 5.5 times after 500 h of water imbibition.

Solute Diffusion Coefficient D e
The value of the solute effective diffusion coefficient can be low due to the high tortuosity of the shale rocks (Barone et al. 1990). Changing the value of D e inherently implies altering the tortuosity of the matrix. It is seen from Fig. 6a that as the value of D e is decreased, damage variable evolves much slower due to hindered solute transport. At 5 h, the damaged area for D e = 3.3 × 10 −16 m 2 /s is significantly lower compared with damaged area for D e = 3.3 × 10 −14 m 2 /s . Similarly, at 50 and 500 h, the damaged area was lower for lower solute diffusion coefficient. Plotting the damage variable evolution at a point in the water-saturated region close to the fracture-matrix interface, (x, y) = (0.02, 7.5) , clearly shows that with decrease in diffusion coefficient, the critical limit of 0.3 is reached much slower leading to delayed micro-fracturing of the rock (Fig. 6b). This implies that micro-fracturing of rock due to chemical damage is delayed significantly in highly tortuous tight rocks.
Since the gradient of the solute concentration between the pore water and fracturing fluid is one of the chemical driving forces causing microstructural alteration of the matrix, the permeability evolution is also significantly slower (Fig. 7). At 500 h, both values of D e lead to a ~ 5.5 times increase in permeability. However, at 50 h, while permeability increased ~ 4.9 times for D e = 3.3 × 10 −14 m 2 /s , it increased ~ 3.9 times for D e = 3.3 × 10 −16 m 2 /s . Similarly, at 5 h, permeability with D e = 3.3 × 10 −14 m 2 /s increased ~ 2.3 times whereas it increased ~ 1.7 times for D e = 3.3 × 10 −16 m 2 /s . Such manifold increase in permeability during the well shut-in periods can cause water entrapment in the micro-structurally altered matrix which cannot be flowed back due to strong wettability forces. These insights highlight the importance of having a multiphysical constitutive theory as ignoring important physical phenomena can obscure actual causes of water uptake in the shale matrix.

Damage Permeability Coefficient ˛D
The damage permeability coefficient D used in Eq. (52) is an experimentally determined parameter. It critically controls the extent of permeability enhancement with damage evolution. Sensitivity of D shows that a slight increase or decrease in value can cause a significant alteration of permeability (Fig. 8). In the preceding discussions, D = 5 was used. When D is increased to 7 (dashed lines in Fig. 8), permeability evolves rapidly with damage with almost an order of magnitude increase in permeability value and a significant increase in the damaged region (here damaged region implies area under the permeability curves). On the other hand, when D is reduced to 2 (dotted lines in Fig. 8), the damaged region is significantly reduced and there is no substantial increase in permeability values. It is noted that the continuum damage mechanics theory is commonly supported by experiments on different materials to characterize the damage variables and a purely theoretical approach is seldom used (Cordebois and Sidoroff 1980;Chow and Wang 1987;Echchur Rangarajan and Ramarathnam 2020). Complex experiments to validate constitutive theories especially for two-phase tight porous systems such as shales remain a challenge to be realized.

Conclusions
We have presented an advanced constitutive framework to account for rock deformation under two-phase flow conditions with time-dependent water-induced chemical damage (microstructural deterioration) to the structure. The Helmholtzian thermodynamic potential was derived with relevant state variables of strain, pore pressures, chemical concentration, and damage. The thermodynamically consistent constitutive framework reproduces Biot's poroelasticity (Biot 1941) and Coussy's unsaturated poroelasticity (Coussy 2010) theories with simple modifications. The constitutive framework based on continuum damage mechanics and the modified mixture theory for multi-phasic porous solids was implemented numerically through solving the problem of shale matrix water uptake using the finite element-based COMSOL Multiphysics platform.
The results show that the theoretical-based first-order damage characterization is more useful than phenomenological or direct approaches (e.g., image-based or measurement of void area ratio) where damage characterization (especially related to micro-fracture initiation) may be inaccurate due to the resolution limit of the imaging equipment. The thermodynamics-based approach adopted here shows that there there is a unique poroelastic interface exists between waterexposed and unexposed regions (David et al. 2015) that was accurately predicted by our damage constitutive theory. One of the main findings was that it is not the magnitude of chemical swelling stresses but the time-dependent chemical damage that causes the microstructural deterioration of the shale matrix (Zhong et al. 2019). This leads to the conclusion that there is substantial permeability enhancement of the shale matrix inducing significant water uptake (and loss). In such a case, it may seem that the entrapment of water in micro-fractures could be the main contributor to fracturing fluid loss into the shale matrix against spontaneous imbibition solely.
The developed constitutive theory can be used to investigate the uptake of the hydraulic fracturing fluid into shale gas formations, design the hydraulic fracturing operation, borehole stability in unsaturated organic-rich shale formations, and also the prediction of gas production and enhanced recovery of shale gas reservoirs. Other practical applications of the theory could be in the context of shallow unsaturated soils, tunnel integrity, nuclear waste disposal, and underground energy storage where variable saturation and chemical damage are expected.

Appendix 1 Balance Laws
The balanced equation for Helmholtz free energy of a thermodynamically open system is given by (Heidug and Wong 1996;Chen and Hicks 2013): where the material time derivative of any argument a is given by: As mentioned in the main text, the energy or thermodynamics-based approach to constitutive modeling is well documented, and hence trivial steps are skipped. For a more rigorous derivation on how to arrive at Eq. (55), the readers are referred to classical literature (c.f. Heidug and Wong 1996;Coussy 2005Coussy , 2010). In Eq. (55), F is the Helmholtz free energy density, is the Cauchy stress tensor, n is the outward unit normal vector, v s is the velocity of the solid, T is the temperature (constant), I w , I g and I c are water, gas, and solute flux, respectively, and is the rate of entropy production density (i.e., the rate of entropy production per unit volume). It is noted that the bold font symbolizes either a vector or a matrix. The flux can be defined as: where = w, g, c for water, gas, and solute, respectively, is the mass density of each component, v is the velocity of each component (water, gas, or solute), and v s is the velocity of the solid phase. The constituents of V in addition to the free energy are the fluid and solid masses which also need to be balanced for any change in time. Since, as already stated earlier, V is closed to solid matter influx and open to fluid (water, gas, and solute) flow, we have: Next, using Reynold's transport theorem and Gauss divergence theorem, we obtain the localized version of the balance equations for free energy, solid mass, and each fluid component mass, respectively: The densities of each fluid component used heretofore were in respect to the unit volume of the bulk (fluid-solid mixture). They are related to their mass densities through: where i=w,g is the porosity of the relevant fluid phase − water or gas. It is noted that the solute component ( = c) is associated with the water phase porosity. The total porosity ( ) of the medium is related to the fluid porosity through the respective fluid saturation (S w,g ): where w + g = 1.

Entropy Production Rate
The entropy production rate in Eq. (55) needs to be specified to quantify the dissipative mechanisms within the rock. It is assumed here that the cause for dissipation is the friction generated during flow at the solid/water interface and due to damage. We assume that the dissipations caused by the inelastic deformation of the rock, and the chemical reactions of the fluid components with one another and with the rock are negligible. In other words, reversible elastic deformation and local interfacial chemical equilibrium are assumed. Hence, using standard arguments of non-equilibrium thermodynamics, we obtain the macroscopic expression for dissipation (Heidug and Wong 1996;Katchalsky and Curran 2014): where D is the damage variable (a scalar), Ḋ its time evolution, and Y D is the thermodynamic force causing damage. The fluxes heretofore were based on the velocity of each component relative to the solid matrix velocity. Following the methodology of Heidug and Wong (1996) and Aghighi et al. (2021), it is more convenient to replace the above fluxes with the diffusion fluxes J that are relative to the mixture's barycentric velocity v f . Using the Gibbs-Duhem equation for multicomponent multiphase systems (Sacchetti 2001) and making necessary rearrangements, the entropy production can be written in terms of two independent fluxes: where u f = S f v f − v s is the Darcy velocity. It is noted that in the above equation, the term ∇( g − w ) can be neglected presuming no interfacial exchange or mixing occurs between the two phases. The readers are referred to Roshan and Aghighi (2012) for details on thermodynamic theory of irreversible processes.

Equations of State
Using the entropy production (Eq. (94)) in Eq. (55) and assuming ∇ ⋅ = 0 (i.e., the rock maintains mechanical equilibrium during deformation), we arrive at: Using the concepts of continuum mechanics and invoking the fluids mass conservation (Eq. 80), we obtain the referential equivalent of Eq. (95) for small strain ( ) problems (Heidug and Wong 1996;Aghighi et al. 2021): Next, assuming that there is no gas residing in the interlayer pore space and compressive pressure is positive, the free energy density of the fluid system in the free pore volume (containing water and gas) can be written as: The above equation can also be written as: To arrive at the appropriate formulation of the skeleton constitutive equations incorporating the skeleton-fluid coupling, it is convenient to adopt a Lagrangian description, i.e., refer the fluid motion to the referential (initial) configuration of the skeleton (Coussy 2005). Invoking the continuum mechanics concepts, we can deduce the referential pore volume free energy density J F pore : (96) F = tr(̇) + wṁw + gṁg + cṁc + Y DḊ . (97) where v w = J w and v g = J g are the pore volume fractions of water and gas per unit referential volume, respectively. The total pore volume is v w + v g . Since J pore = m pore , we have: Subtracting this from the combined shale/fluid bulk system, the free energy of the 'wetted' rock skeleton (i.e., solid part of the rock including the clay interlayer water) is obtained (c.f. (Ma et al. 2022) for more details): Here, m bound = m − J pore is the mass density of the = w, c fluid component in the interlayer pores in the reference volume. It is reiterated that there is no such 'bound' gas in the interlayer pores or the 'wetted' matrix which is thermodynamically consistent with our earlier assumption that gas adsorption is trivial (or negligible at least in inorganic constitutes that form majority of the matrix volume). Next, it is convenient to employ the following potential based on Eq. (101): Using the local state postulate, the rock skeleton thermodynamic potential can be written in the form (Coussy 2005): The derivative of W , as a function of the state variables (after appropriate Legendre transforms in v w , v g , and m bound , can then be written as:

Appendix 2 Damage Evolution
As discussed in the main text, when microstructural deterioration of material is discussed in the framework of continuum mechanics, the internal state variable called the damage variable (D) , is utilized. In this study, this variable is considered isotropic across the REV and hence is a scalar. (99) (103) W = W , D, p w , p g , .
The microstructural deterioration will cause only slight micro-cracking which does not necessarily lead to permanent deformation. However, such deterioration is manifested in the reduction of mechanical properties of the material due to a reduction in grain-to-grain cohesion. In other words, the energy dissipated due to microstructural deterioration is used to reduce the material cohesion by chemical reaction influencing the mechanical properties although a minor initiation of the micro-cracks due to such chemical reaction is inevitable. The strain regime can remain elastic during such deterioration of the material which is true for most brittle materials (Khan et al. 2004(Khan et al. , 2007Murakami 2012;Jia et al. 2020). Also, since there is no cyclic mechanical loading nor chemical drying-wetting cycles, this consideration of negligible permanent deformation remains valid. The dissipated energy due to microstructural deterioration can be represented as the loss in the elastic strain energy density w e of the REV. The elastic strain energy density is defined in the stress-strain space as (Lemaitre 1992): It is recalled that the thermodynamic force causing the microstructural alteration recalled from Eq. (15) is Y D . Based on the Helmholtz free energy potential (W) , it is interpreted as (Lemaitre 1992): which physically implies that the measure of skeletal free energy change due to any change in microstructural deterioration (damage). It is noted that, as per the choice of definition, Y D can be equated to a quantity −Y D to work with a positive quantity (Lemaitre 1992) (the readers are advised that Y D is denoted by Y in Lemaitre (1992)). The thermodynamic force of damage Y D can be related to the elastic strain energy density by (Lemaitre 1992): The above equation implies that the thermodynamic force associated with damage Y D is proportional to the elastic strain energy density with a proportionality constant 1∕(1 − D) . The proportionality constant increases with an increase in damage (D) . Thus, if the strain energy density w e can be calculated, the thermodynamic force causing damage can also be subsequently estimated by substituting w e in Eq. (107). In the stress-strain space, Eq. (105) becomes: (108) For chemical damage due to the combined action of water-solute-clay interactions, we use the similar concept of strain energy density to nurture the expression for the evolution of D as discussed in the main text.

Model Validation
The full set of coupled field equations cannot be solved together analytically due to high nonlinearity. Hence, decoupled physics (i.e. a sub-model) is separately solved and presented to demonstrate the accuracy and robustness of the developed constitutive theory. Validation of decoupled physical parts of the system is provided by comparison with either analytical or previous literature data.

Validation of Poroelasticity
The poroelasticity physics is validated through the reproduction of the example by Zheng et al. (1955) and Freeman et al. (2008) for a rock sample undergoing uniaxial compression (i.e., a drained uniaxial single-phase poroelastic case is solved). This example consists of a 2 m × 3 m 2D porous medium fully saturated with oil where only the top boundary is open for fluid to flow (Fig. 9). A 4 MPa uniform load is applied on the top boundary while the other sides have roller constraints. The input data used are shown in Table 2 2. To solve for simple poroelasticity, Eqs. (20) and (41) are utilized without the gas phase, and solute concentration terms, and Eqs. (42), (46), and (51) are made redundant.
The transient pore pressure and vertical displacement profiles are shown in Fig. 10 for different times. The results are an excellent match with the analytical and numerical results presented by Zheng et al. (1955) and Freeman et al. (2008). Furthermore, the effective stress is plotted in Fig 10Fig. 11 which shows that with time, the effective stress approaches the applied load as the pore pressure dissipates.

Validation of Two-Phase Flow
The analytical solution for counter-current spontaneous imbibition is presented and compared with the solution given by COMSOL. The boundary conditions, input parameters, and simulation geometry and domain are the same as those presented in Sect. 3.2. It is noted that the analytical solution is limited to the infinite-acting domain, i.e., before water saturation reaches the boundary at the other end of the REV. The analytical solution is 'semi-analytical' in the sense that it is an integral solution requiring trivial iterations to Since a given saturation value travels a distance proportional to √ t , the water flux scales as 1∕ √ t with a proportionality constant C which has units of m∕ √ s. The magnitude of C determines the speed of spontaneous imbibition and hence is analogous to the diffusion capillary coefficient ( D c ). The solution is arrived at by introducing a new fractional-flow function ∀ (Schmid et al. 2011) such that: ∀ here is distinct from the Buckley-Leverett fractional flow, although is inspired from it. Here, the fractional-flow function is represented by ∀ to distinguish it from the Helmholtz free energy symbol F . The final equation which is solved is a second-order ordinary differential equation for ∀: Apart from porosity, permeability, water/gas viscosities, the required inputs are S wi , S or , k rw , k rg , and P c which are provided through fractal theory. The water saturation profiles obtained from the simulation merge into one curve when scaled by x∕ √ t . This scaled curve is then matched with the water saturation versus similarity variable ( ) curve obtained from the analytical solution. It can be seen (Fig.12) 12) that an excellent match is obtained between simulated and analytical profiles. The value of C obtained after iterations was 6.03x10 −9 m∕ √ s . This validation endorses the accuracy of the two-phase flow formulation of the constitutive theory.

Validation of Solute Transport
The validation of the solute transport equation is provided by comparing the solute mass fraction profiles from COM-SOL with those from Heidug and Wong (1996). To solve for solute transport only, Eq. (46) from the set of constitutive equations is considered, and Eqs. (20), (41), (42), and (51) are made redundant. For the solution, the geometry and boundary conditions enacted by Heidug and Wong (1996) are used (Fig. 13). The geometry consists of a 2D rectangular domain where the left side of the saturated sample is exposed to a brine containing 0.35 mass fraction ( C c x=0 ) of salt. The inside of the sample is initially saturated with distilled water without any salt, i.e., C c t=0 = 0 (or C inside = 0 (111) ∀∀ �� = − � 2C D c .   The results are plotted in Fig. 14 which are along any line parallel to the X-direction. The solute mass fraction profiles obtained from COMSOL simulation are represented by solid lines whereas those from Heidug and Wong (1996) are represented by dashed lines. It can be seen that solute mass fraction profiles obtained from COMSOL simulation at different times are very close to those reported by Heidug and Wong (1996). The apparent differences in the curves can be attributed to the linearization of the relationship between solute chemical potential ( c ) and solute concentration (C c ) in this study. On the other hand, the relation between c and C c . was logarithmic in the study by Heidug and Wong (1996). Nonetheless, as can be seen, the differences are minimal.

Fig. 13
Geometry and boundary conditions for solute transport from Heidug and Wong (1996)  in Fig. 13). It is noted that C c here represents the solute mass fraction which is unitless. Other parameters used in the simulation are tabulated in Table 33.