A hydrodynamic model for chemo-mechanics of poroelastic materials

Geomaterials are susceptible to dissolution, whereby solid particles disintegrate into surrounding aqueous fluids. The dissolution rate depends on multiple factors, such as grain mineralogy and morphology, hydraulic and thermal environments, and stress field (Gautelier et al., 1999; Xie et al., 2011; Thyagaraj & Das, 2017; Sun et al., 2018). As the solid skeleton of geomaterials dissolves, their strength weakens and stiffness softens (Shin & Santamarina, 2009; Stefanou & Sulem, 2014; Shen et al., 2019; Viswanath & Das, 2020). Based on mineralogical composition, geomaterials can either dissolve into water through physical processes such as salt dissolution (Shalev et al., 2006) or react in acidic solutions through chemical processes such as the dissolution of carbonate rocks or calcareous sands (Kalia & Balakotaiah, 2007). Owing to the widespread distribution of carbonate minerals in nature (Witze, 2013), solid dissolution induced by chemical reactions has attracted extensive scientific attention in the geotechnical and geological communities (Morse & Arvidson, 2002). In practice, chemical dissolution through acid injection can be leveraged as a stimulation technique to enhance the productivity of wells drilled into carbonate rocks; this technique is exploited in the oil and gas industry and geothermal reservoirs (Schumacher & Schulz, 2013; Hu & Hueckel, 2019). However, it can also undermine the stability and sustainability of geotechnical and geological systems, with the opening of sinkholes (Basso et al., 2013), breach in carbon dioxide (CO₂) sequestration (Dawson et al., 2015) and the triggering of landslides by acid rain (Xiao et al., 2021).

A number of approaches have been proposed in the past to explore and quantify the impact of chemical reactions on mechanical constitutive features. In one approach, the well-known Cam-Clay model was extended by incorporating a new weathering index to control plastic hardening and approximate the effect of chemical dissolution on mechanical properties (Nova et al., 2003; Buscarnera & Das, 2016). A different approach examined the connections among various material features through a three-scale model (Hu & Hueckel, 2007; Ciantia & Hueckel, 2013), where the relative mass removal of dissolvable solids was used as a measure of weathering. This approach was further adjusted to analyse the stiffness degradation of calcarenite samples both before and after the elastic limit (Gajo et al., 2015; Ciantia & Di Prisco, 2016). These models are able to consider the intrinsic connections between the micro and macro scales with the objective of reproducing the complex responses of calcarenite samples subjected to various testing conditions.

Although the previously adopted elastoplastic theory has provided useful predictions, corresponding evolution laws have not always been grounded on rigorous thermodynamics, including the first and second laws, Onsager's reciprocity restrictions and material objectivity. Proper considerations of thermodynamics have the benefit of eliminating over-fitting and keeping model constants to a minimum. Another benefit of thermodynamic-based models is that they can be used to study previously untested experimental conditions with confidence. A most complete and rigorous way to account for all the fundamental subtleties implied by thermodynamics (beyond just first and second laws) is to derive constitutive equations using the hydrodynamic procedure. However, this comprehensive approach has not yet been developed to enable modelling of the chemo-mechanically coupled behaviour of poroelastic geomaterials.

The aim of the current paper is to fill this theoretical gap. Thanks to the rigorous derivation, the resulting model will be shown to require only physically motivated and clearly stated assumptions, associated with only meaningful and much fewer parameters than have ever been used before for these problems. The derivation follows the hydrodynamic procedure, which was pioneered by Landau (Landau & Lifshitz, 1980, 1987), and extended to address super fluidity (Khalatnikov, 2018) and liquid crystals (De Gennes & Prost, 1993). Recently, Jiang & Liu (2009) adopted it for granular solid systems, with a number of later developments, all focusing on sand (Zhang & Cheng, 2016; Alaei et al., 2021; Zhang & Soga, 2021). This theory was further developed for fully and partially saturated soils under both equilibrium and non-equilibrium conditions (Jiang et al., 2017; Einav & Liu, 2018, 2022). Models developed from the hydrodynamic procedure cogently satisfy the thermodynamics and conservation laws as well as symmetry and compatibility restrictions (Jiang et al., 2017; Einav & Liu, 2018), thus avoiding many possible physical inconsistencies. To hydrodynamically accommodate solid dissolution, a key step in the current derivation is to incorporate mass fractions of all the ionic species in the pore fluid as internal state variables and capture the stiffness degradation in terms of density changes. In this paper, the considered model exhibits a linear poroelastic response only in the absence of chemical reactions. However, upon chemical reactions, that same material responds non-linearly due to density and thus stiffness changes.

Aiming to be both general and specific, this paper is organised as follows. The theory is developed for any solid–fluid mass exchange process. The physical and mathematical grounds thereby established lay down the general seeds for the comprehensive modelling of the role of solid dissolution on material deterioration. For the sake of simplicity, the model specified later considers a poroelastic solid medium and focuses on understanding the influence of chemical dissolution on the elastic properties of calcium carbonate under fully saturated conditions. This theoretical simplification can be directly applied to study chemo-mechanical effects in porous rocks prior to their elastic limit. In this context, the proposed model is adopted to interpret tests on calcarenite subjected to acid invasion as undertaken by Ciantia and colleagues (Ciantia et al., 2013, 2015; Ciantia & Hueckel, 2013). Complex features such as inelasticity at higher pressures and deformation-dependent chemical kinetics are left for future considerations. Even without those effects in the current model, the pure elastic formulation of this paper can still successfully explain a number of observed phenomena using much fewer parameters than ever before, all being physically meaningful.

In order to help readers see the complete picture before entering into specific details, a brief road map is first described in Fig. 1. In particular, this road map helps to visualise the essential steps in the current hydrodynamic modelling of geomaterials subjected to a chemo-hydro-mechanical influence. In that diagram general hydrodynamic concepts as well as more specific assumptions for saturated poroelastic materials are denoted by regular and bold fonts, respectively. It is noted that while the integration of all the general hydrodynamic concepts is indeed rather delicate mathematically, it follows a standardised physical procedure, which only needs to be done once. Once established, as shown later, the rigorously derived relationships can be used directly to capture the physico-chemical behaviours of interest and model the dissolution processes of geomaterials, including the material non-linearity, temperature dependence and the coupling of the various transport processes. Future modelling effort should mainly focus on refining the conserved energy function and Onsager's transport matrix based on the material features of interest, since all the other constitutive relationships can be systematically recovered from the consolidated relationships in this paper. Although some additional mathematics would be used to interpret key mechanisms, they are directly recovered from the stated energy potential based on well-established physical observations, instead of treating them inconsistently in a disconnected phenomenological way, as done in previous presentations.

As is shown, the hydrodynamic procedure always starts with the definition of the conserved energy (U). It is a function of a series of independent state variables representing the system properties and their energy-conjugate counterparts. Each state variable follows its own balance law, such as the momentum balance, entropy balance and mass balance. Combining these with energy conservation (i.e. first law of thermodynamics), one can derive the entropy production, which should always be non-negative to preserve the second law of thermodynamics. To guarantee this, the Onsager's reciprocal theory (Onsager, 1931) is typically leveraged to relate ‘velocities’ to ‘forces’. Consequently, considering loading and boundary conditions, the field equations characterising the chemo-hydro-mechanical coupling can be obtained by integrating the balance equations with the constitutive relations (i.e. stress–strain and Onsager's relationships) and used to analyse the systematic behaviour of geomaterials.

This section lays out the general structure of the hydrodynamic procedure adopted for geomaterials subjected to a chemically active environment.

As part of the hydrodynamic procedure (Landau & Lifshitz, 1980, 1987), the conserved energy density is first defined in both the moving and rest (i.e. zero barycentric velocity v_i = 0) frames as U and u, respectively. These energies are linked through the kinetic energy k

where u is also known as the internal energy density; g_i = ρv_i represents the momentum density; ρ is the partial density of the continuum; and i is the Cartesian coordinate, with implicit summation of repeated indices.

Theoretically, the conserved energy density potentially depends on many independent internal variables that most generally describe the state of the material. However, to achieve clarity it is crucial to start by scoping down the description using a sufficient set of simple internal state variables, with which the complex state of the material could be represented while preserving the capability to capture key features of interest. As discussed by Einav & Liu (2018), the conserved energy density for porous materials could be generally considered as

where s represents the entropy; ε^e_ij is the elastic strain tensor; and ρ_β is the partial density of the involved β phases. These partial densities are specified (Fig. 2) for the solid (s) and the aqueous fluid (f) for a fully saturated system (i.e. β ∈ {s, f}), or for the solid (s), aqueous fluid (f) and air (a) for a partially saturated system (i.e. β ∈ {s, f, a}), or for the solid (s) and the air (a) for a fully dry system (i.e. β ∈ {s, a}). Within a representative volume (V) with a total mass (m), each phase has its own mass (m_β), volume (V_β), and partial (ρ_β) and intrinsic (⁠$ρ^β$⁠) densities, so that

from which the β-phase volume fractions (ϕ_β) and their sum can be defined

Owing to negligible ionic concentrations in natural water, the pore fluid is generally considered as pure water (f ∈ {w}). However, the occurrence of solid dissolution can significantly change the ionic concentrations in the pore fluid, or even introduce new ionic species. Therefore, instead of pure water, the pore fluid is an aqueous solution composed of water and all the ionic species (f ∈ {w, γ}) (Fig. 2)

where ρ_w and ρ_γ are the partial densities of water and γ-ionic species defined in terms of the total volume V (i.e. ρ_γ = m_γ/V, ρ_w = m_w/V, with m_w and m_γ being their corresponding masses), with $ρ~w$ and $ρ~γ$ being the partial densities in the solution characterised in terms of fluid volume V_f (i.e. $ρ~γ=mγ/Vf$ and $ρ~w=mw/Vf$⁠). Therefore, the mass fractions of each species in the pore fluid can be defined as

where c_w and c_γ are the mass fractions of the water and the ionic species with respect to the mass of pore fluid, respectively. Thus, the molar concentration of the ionic species x_γ could be computed as

where M_γ represents the molecular mass of the γ species.

Each component in the pore fluid requires one independent internal variable to describe its presence. Theoretically, it can be either mass fraction, density, or concentration. Here ρ_f and c_γ are selected as independent variables while all the other quantifies can be analytically determined from them. Accordingly, the conserved energy and internal energy densities are given as

Following the hydrodynamic procedure, the energy-conjugate variables for each of the internal state variables could be defined as (Jiang & Liu, 2009; Jiang et al., 2017; Einav & Liu, 2018)

where T is the thermal temperature; μ is the overall chemical potential; μ_β is the partial chemical potential of the β-phases; and X_γw is the chemical potential difference between the γ-ionic species (μ_γ) and water (μ_w) defined as

while σ^e_ij is the elastic stress defined as (Einav & Liu, 2018; Alaei et al., 2021)

with δ_ij being the Kronecker delta, σ_ij and σ^D_ij the total and viscous stress tensors, respectively, and p^T the thermodynamic pressure, a well-known property in physics that is given by the rate of energy with respect to volume (i.e. p^T ≡ ∂(UV)/∂V). For a unit mass of soil mixture (i.e. V = 1/ρ), the thermodynamic pressure is computed as (Jiang & Liu, 2009)

The sum of momentum densities of all the components is equal to the total momentum density of the mixture

where v^β_i, v^w_i, v^f_i and v^γ_i are the flow velocities of the corresponding components. The total chemical potentials of the entire mixture (μ) and the aqueous fluid (μ_f) could be eliminated from equation (9), as well as equations (3) and (5)

Each of the previously defined state variables requires its own evolution law. Specifically, the momentum density is described by its conservation law

where G_i is the gravitational acceleration.

The entropy has its own balance law, while the rate of entropy production $R$ is always non-negative by virtue of the second law of thermodynamics

where sv_i and f_i are the convective and dissipative entropy currents, respectively.

The elastic strain tensor is obtained from (Einav & Liu, 2018)

where $ε˙ij=−(∇ivj+∇jvi)/2$ and $ε˙ijp$ are the total and plastic strain rate tensors, respectively; v_k∇_kε_ij^e represents the advection of elastic strains; while $Ωij=(∇ivj−∇jvi)/2$ is the anti-symmetric component of the velocity gradient. Notice that the term $εikeΩkj−Ωikεkje$ accounts for the rigid body rotation of the elastic strain, which satisfies

Finally, the evolution of partial densities and mass fractions of the ionic species are given by their mass balances. Appendix 1 provides the general formulation for systems where the overall mass is being conserved, while mass transfers due to multiple parallel chemical reactions are being balanced for the individual components. Here, for simplicity, the result is specified for systems undergoing a single chemical reaction

where M_βN_βξ and M_γN_γξ act as source/sink terms for the density evolution of the β phases and γ species attributed to the chemical reaction; M_β and M_γ indicate the molecular mass of the β phases and γ species; ξ (with units of mol/(m³h)) quantifies the rate of chemical reaction governed by chemical kinetics; N_β and N_γ represent the stoichiometric numbers of the β phases and γ species, which are equal to zero for components not involved in the reaction. Conventionally, the stoichiometric numbers take negative values for reactants and positive values for products. Finally, J^β_i and J^γ_i are the dissipative density fluxes of the involved components (β or γ) relative to the barycentric velocity, which are being defined as

Combining equations (3), (5), (13) and (20)

Based on the hydrodynamic procedure (as detailed in Appendix 2), the non-negative rate of entropy production could be derived

where $A=μp+μr$ represents the sum of the reactants’ (μ_r) and products’ (μ_p) chemical potentials, known as the driving force of the chemical reaction, recalling that the signs of μ_p and μ_r are opposite due to the sign convention of the stoichiometric numbers. The resulting sign of $A$ determines the direction of the chemical reaction, with negative values signifying forward reaction (from reactants to products), positive implying reverse reaction and zero representing chemical equilibrium (Job & Rüffler, 2016). In the absence of chemical reaction (i.e. J^γ_i = 0 and $A=0$⁠) the derived $R$ naturally converges to the mathematical result presented in Einav & Liu (2018). Here, the occurrence of chemical reaction contributes to the entropy production through additional terms from chemical kinetics (⁠$−ξA$⁠) and the dispersion of ionic species (⁠$−∑Jiγ∇iμγ$⁠) (Gao et al., 2021).

By potentially considering the general scenarios of solid dissolution and hydromechanical conditions (i.e. variable degrees of saturation and different inelastic processes, such as grain crushing and rearrangement), the proposed modelling framework can be specified to analyse any dissolution process in both saturated and unsaturated soil systems with the flexible incorporation of advanced constitutive models. This theoretical generality offers a way to quantify the role of chemo-mechanical coupling during the solid dissolution in any complex geomaterial and investigate the underlying physics of material response, which could benefit the assessment, prevention and mitigation of natural disasters associated to solid dissolution.

The derivation in the previous section was rather general. Conversely, the following discussion is limited to the study of solid dissolution in saturated poroelastic materials composed of an incompressible solid phase (i.e. β ∈ {s, f}, $ε˙ijp=0$ and $ε˙ij=ε˙ije$⁠). More specifically, the primary focus is on the dissolution of calcium carbonate in an acidic environment prior to its elastic limit, at isothermal room temperature (i.e. T = 24°C, $∇iT=0$ and $T˙=0$⁠), and under quasi-static conditions of stationary samples in the laboratory (i.e. negligible v_i, U = u and σ^D_ij = 0). Such simplifications exclude any physical intricacy that may arise from material inelasticity and surface tensions, while highlighting the critical role of the chemo-mechanical processes in the elastic characteristics.

The internal energy density reflects the total energy preserved by all the phases in the geomaterial continuum, which can be calculated as the sum of the elastic strain energy density preserved by the microstructural skeleton u^e, and the free energy density of the solid–fluid mixture u^m possessed by the involved single-phase domains, which controls the thermodynamic pressure

In particular, the elastic strain energy is formulated to produce linear elasticity, with density-dependent scaling of stiffness motivated by (Gibson & Ashby, 1982; Rubin & Einav, 2011; Alaei et al., 2021)

where ρ_s* is the unstressed solid density, which in general geotechnical settings remains approximately equal to the intrinsic counterpart of the variable (i.e. $ρs∗≈ρ^s$⁠, therefore ρ_s/ρ_s* ≈ ϕ_s, with ϕ_s being the solid fraction); $K¯$ and $G¯$ are the bulk and shear stiffness constants satisfying $G¯=3K¯1−2ν/21+ν$⁠; ν is the Poisson's ratio; and n is a power coefficient reflecting the density-dependent elasticity, which for foams and soils is typically taken as 3 (Gibson & Ashby, 1982; Viggiani & Atkinson, 1995). Rather than tuning this coefficient against experiments, this parameter is here fixed to 3, to ensure predictions remain genuine, rather than acting as curve-fits, as done in the past. Thus, by taking n = 3 and by using equations (9) and (23), the corresponding elastic triaxial stresses become

where p^e is the mean elastic stress and q^e is the deviatoric elastic stress, which here, in the absence of a viscous stress equals to the deviatoric stress (q = q^e). Upon further differentiation, the instantaneous elastic bulk (K) and shear moduli (G) are obtained

so as the solid skeleton dissolves and the solid fraction decreases, the elastic stiffness may soften non-linearly, albeit having the elastic stresses (equation (25)) linearly depend on the elastic strains.

The free energy accounts for the volumetric compressibility and osmotic concentration of the two large single-phase domains made of solid and fluid, which could be calculated analytically (Jiang et al., 2017; Einav & Liu, 2022)

where $u^s$ and $u^f$ are the intrinsic free energy densities of the individual phases of the solid and fluid, respectively. It is worth noting that for incompressible solids (i.e. K ≪ K_s), the intrinsic solid density (⁠$ρ^s$⁠) is nearly constant for typical geotechnical problems, in which case the intrinsic fluid density (⁠$ρ^f$⁠) directly depends on the partial densities based on equation (4). Therefore, the intrinsic densities in equation (27) are not regarded as independent state variables. However, subtle changes to the intrinsic solid density $ρ^s$ still need to be accounted for in order to recover the intrinsic thermodynamic pressures and intrinsic chemical potentials (Jiang et al., 2017), respectively

For the fluid domain the intrinsic thermodynamic pressure physically represents the total energy that constitutes the pressure head, which is equal to the hydrostatic pressure used in mechanics. In other words, $p^fT$ can be seen as the hydrostatic pressure of the pore fluid and $p^sT$ as the internal pressure within the solid constituting the particles as they are being compressed isotropically by the surrounding fluid. For saturated cases, since the solid domain is solely compressed by aqueous fluid, $p^sT=p^fT$⁠, which is conveniently marked as the common pressure hereafter (i.e. $P=p^sT=p^fT$⁠). Such equality can also be theoretically proven through energy minimisation as done by Jiang et al. (2017). This proof, however, involves additional mathematical derivation, which is skipped to avoid repetitions and further complications to the current paper. Curious readers may follow this derivation in Jiang et al. (2017).

As presently stated in equation (12), the thermodynamic pressure of the mixture p^T depends on both u^e and u^m through u. However, since the elastic energy density scales as u^e ∝ (ε^e_ij)², and the elastic strains are normally small with the order O(ε^e_ij) ≪ 1, the contribution of u^e to p^T is negligible. In contrast, the contribution of u^m, which is the dominating term controlling the thermodynamic pressure, can be of any magnitude (Alaei et al., 2021; Einav & Liu, 2022) so that (Appendix 3)

Considering equations (11) and (29), the elastic stress for stationary systems (i.e. σ^D_ij = 0) could be identified as Terzaghi's effective stress

which is valid for the case of incompressible solid constituents given the current choice of the internal energy function in equation (23). For compressible solids (Lade & De Boer, 1997; Gajo, 2010), this energy needs to be generalised to accommodate the changes in the intrinsic free energy (u^m) and the elastic strain energy (u^e), keeping in mind the Biot's effective stress principle (Biot, 1941; Nur & Byerlee, 1971; Detournay & Cheng, 1993)

where K_s is the intrinsic bulk modulus of a pure solid typically ranging from 20 to 70 GPa (Madhubabu et al., 2016; Sohn et al., 2017). However, such a generalisation is outside the scope of the current derivation, since for the soft calcarenite studied in this work, solid compressibility is relatively negligible owing to its relatively high porosity (K ≪ K_s), as demonstrated throughout the model evaluation. Note that in this limit of K ≪ K_s, the Biot's relation reduces back to Terzaghi's relation.

To quantitatively evaluate p^T, it is essential to formulate the intrinsic free energy analytically for each single-phase domain. Specifically, the free energy of the solid part is formulated to capture the primary features of the solid domain: (a) its volumetric compressibility, and (b) its pressure-dependent chemical potential. Considering these two conditions, this free energy is taken as

where the superscript ⁰ indicates variables at the standard state (i.e. room temperature T = 24°C, atmospheric pressure P⁰ = 101 kPa and a molar concentration of x⁰_γ = 1 mol/l for aqueous solution). Based on equation (28)

which converges to the standard values at the standard state. Moreover, the intrinsic density of a solid varies very slightly in geotechnical practices (i.e. $ρ^s/ρs0→1$⁠). By ignoring the higher order terms of the Taylor series

therefore, the exact chemical potential of the solid in equation (33) can be approximately expressed as

which is a widely recognised relationship in physico-chemistry to assess the pressure dependency of the chemical potential (Job & Rüffler, 2016). Taking the increments of equation (33)

which recovers the definition of solid compressibility (Jiang et al., 2017). In other words, the volumetric compressibility is the dominating mechanism provoking the changes in chemical potential and pressure as designated when formulating the intrinsic free energy of the solid domain.

For pore fluid the intrinsic free energy needs to incorporate the osmotic concentration due to the presence of ionic species, in addition to the volumetric compressibility. Here, the pore fluid is assumed to be a substantially dilute ionic aqueous solution (x_γ ≪ 1 mol/l). The volumetric compressibility is mainly attributed to the solvent, while the solute primarily impacts the osmotic concentration (Elliott, 2020). In this context, the intrinsic energy is separated into two terms to restore the two mechanisms

where K_w = 2 × 10⁶ kPa is the intrinsic bulk modulus of water and R = 8·314 J/(K mol) is the ideal gas constant. It is important to note that although equation (37) is expressed in terms of $ρ~w$ and x_γ, they are not the independent state variables chosen to compute the internal energy, but are functions of the selected state variables $ρ^f$ and c_γ as shown in equations (5) and (7). Based on equation (28), the pore fluid pressure becomes

with its increment being

in which the first term on the right-hand side of equation (38) reflects the pressure contribution of fluid compressibility (⁠$P~w$⁠) provided by the solvent (i.e. water), while the second term characterises the osmotic pressure (⁠$∑P~γ$⁠) induced by the dissolved ions, which converges to the Van't Hoff equation (Thyagaraj & Salini, 2015). In other words, the pore fluid pressure is simultaneously governed by the competing mechanisms of water compressibility and osmotic effect. For example, imagine adding salt into a cup of water to create seawater with a salinity equivalent to 35 g salt per litre of water, as shown in Fig. 3. Since the fluid pressure should maintain a constant atmospheric pressure throughout the dissolution (i.e. P = P_a), the pressure variation due to compressibility and osmotic effects should cancel out (i.e. $dP~w=−∑dP~γ$⁠). Based on equation (39), the calculated seawater volume (V_f) is 1·0015 litres, which is higher than the initial water volume (V_w = 1 litre) but lower than the sum of water volume and salt volume (m_s = 35 g, $ρ^s=2160$ g/l, V_s = 0·0162 litres, and V_w + V_s = 1·0162 litres). Then, the intrinsic density of seawater is computed as $ρ^f=1033$ g/l, which is the typically reported density of seawater at the room temperature. This suggests that the generated osmotic pressure from the increase of ionic concentrations during salt dissolution is cancelled by the pressure reduction due to water compressibility as the fluid volume increases and, consequently, the water density in the solution decreases (⁠$ρ~w=998⋅5$ g/l).

Furthermore, based on equation (9)

Therefore, according to equation (10) the chemical potentials of the solute (i.e. ionic species) and solvent (i.e. water) are

As the water density in the solution is also close to its quantify at the standard state (i.e. $ρ~w/ρw0→1$⁠), the use of a similar approximation to that in equation (34), as well as equations (38) and (41), suggests that the exact chemical potential of the water in the solution can be approximated by

which is a widely recognised relationship in physico-chemistry to account for the chemical potentials of both solute and solvent in a given solution (Atkins & De Paula, 2014; Job & Rüffler, 2016; Elliott, 2020).

As stated at the beginning of this section, the studied system is composed of saturated poroelastic material under isothermal and quasi-static conditions (i.e. $∇iT=0$⁠, β ∈ {s, f}, σ^D_ij = 0, and $ε˙ijp=0$⁠). In this context, considering equations (10) and (21), the rate of entropy production in equation (22) reduces to

where X_ws = μ_w − μ_s quantifies the chemical potential difference between water and solid. Considering the second law of thermodynamics, the entropy production $R$ should always be non-negative. To satisfy this constraint, most generally one may adopt the full Onsager's reciprocal theory to allow coupling among all the dissipative terms involved above for microscopic processes (Onsager, 1931). However, the determination of the coupled Onsager's transport coefficients requires delicate experiments so that an adequate constitutive calibration could be carried out. Therefore, in the absence of such experiments, here, for simplicity, a linear proportionality law is adopted to describe the relationships between the three ‘forces’ involved (as termed by Onsager, i.e. $−∇iXws$⁠, $−∇iXγw$ and $−A$⁠) and the corresponding ‘velocities’ (as termed by Onsager, i.e. J^f_i, J^γ_i and ξ)

which readily ensures the positiveness of the entropy production in equation (43) given positive proportionality coefficients. As proved below, such proportionality relations can automatically recover the well-known empirical laws of Darcy's law, Fick's law and the law of chemical kinetics.

The proportionality between J^f_i and $∇iXws$ is expressed as

where r_f is taken as a positive linear coefficient for fluid diffusion. Therefore, considering equations (29), (39) and (41)

For systems where the pore fluid is solely water (i.e. x_γ = 0, ρ_f = ρ_w, $ρ~w=ρ^w$ and $Jif=Jiw$⁠), the last term representing the osmotic effect disappears as $u~γ$ vanishes and only water compressibility persists, therefore

while using equations (13) and (20)

which converges to Darcy's law with k_w being the hydraulic conductivity of the water.

The proportionality between J^γ_i and $∇iXγw$ is expressed as

where r_γ is the positive diffusion coefficient of the γ-ionic species. This law suggests that apart from diffusion induced by pressure gradients, the spatial gradient of ionic concentrations also promotes diffusion. Combining the above with equation (40), the ionic diffusion driven by the concentration gradient can be written as

Under the assumption of incompressible fluid (i.e. $∇iρ~w=0$⁠) the above relation reduces to provide the famous Fick's law of diffusion

where d_γ is the diffusion coefficient.

The proportionality between ξ and $A$ is expressed as

where r_c is a positive state-dependent reaction coefficient controlling the chemical kinetics of the chemical process of interest. For example, this work focuses on the chemical dissolution of calcium carbonate (CaCO₃), which is typically present in calcites, calcarenites, limestones and calcareous sands. This solid mineral is here considered to be fully saturated in an acidic solution (H⁺), and together these are the system's reactants. With sufficient presence of hydrogen ions (H⁺) in the pore fluid (typically, when the pH is less than 4), the overall reaction kinetics is dominated by the forward reaction (i.e. $A<0$⁠) below (Hoefner & Fogler, 1988; Li et al., 2008; Boyd, 2019)

where the system's products are given by Ca²⁺ as calcium ions, H₂O as water and CO₂ as carbon dioxide. Note that the backward reaction (i.e. $A>0$⁠) is also a naturally important process, which represents mineral precipitation (Prigiobbe et al., 2009). Although nothing prevents the theory from capturing mineral precipitation, this is outside the scope of the current work. Furthermore, it is assumed that the carbon dioxide produced can escape the system immediately after generation, without evoking surface tension (Fernandez-Merodo et al., 2007), but it does play a role in influencing the mass conservation in equation (19) and as a driving force of the chemical reaction in equation (22). Consequently, the pore space is fully saturated with an aqueous solution composed of dissolved calcium salt and acid. Although various types of anions and cations might arise in the pore fluid, such as Ca²⁺, H⁺, HCO₃⁻ and OH⁻, only those directly taking part in the reaction described in equation (I) are considered in the current framework (γ ∈ {H⁺, Ca²⁺}), while the changes in the concentrations of other ions are relatively negligible. In this context, the molecular mass, stoichiometric number and chemical potentials at the standard state of the species involved are listed in Table 1. As the pressure and ionic concentrations shift away from the standard state, the value of chemical potential needs to be updated correspondingly based on equations (35) and (41), which, consequently, leads to an evolving $A$ and indirectly introduces a source of pressure and concentration dependency to the reaction kinetics.

Furthermore, it is widely understood that in strongly acidic environments (pH lower than 4) the kinetics of calcite dissolution is primarily controlled by the transporting mechanism of hydrogen ions in and out of the reaction surface (Berner & Morse, 1974; Morse & Arvidson, 2002). Therefore, any property that can perturb pore fluid diffusion would eventually modulate the dissolution kinetics of carbonate minerals. This includes features associated with pore morphology (i.e. pore volume, surface area, heterogeneity, etc.) and fluid properties (i.e. flow rate, temperature, pH, etc.) (Berner & Morse, 1974; Morse & Arvidson, 2002; An et al., 2021). To capture the experimental findings of the calcite dissolution kinetics, various geochemical models have been proposed to address the impact of different dependent features, among which the reaction rate is normally a function of the molar concentration of hydrogen ion x_H (i.e. pH value) of the bulk fluid under isothermal conditions (Gautelier et al., 1999; Li et al., 2006; Ciantia & Di Prisco, 2016). Therefore, the reaction rate is here quantified by the following state-dependent reaction coefficient

where b is a power coefficient and k_c > 0 (with units of kg mol/(J m³ h)) is the reaction rate constant. Experimental studies suggest that the rate of calcite dissolution is linearly proportional to the hydrogen ion concentration at low pH regimes (i.e. b ≈ 1 for pH < 4) (Fredd & Fogler, 1998; Morse & Arvidson, 2002), but the problem becomes more complex when the pore fluid transits from the acidic to the alkaline regime (pH > 4), in which case the reaction kinetics is governed by both mass transport and surface reaction (Brantley et al., 2008). Therefore, here b is treated more generally as a phenomenological model parameter to approximate the reaction kinetics in both strong and weak acid environments. Furthermore, the model lumps together all the possible mechanical effects on the chemical kinetics into the calibration of the single reaction rate coefficient k_c, despite the fact that physically, this coefficient is sensitive to pore structure since the change of surface area can significantly alter the reaction kinetics. However, obtaining such a state-dependent reaction kinetics would require an expensive and careful experimental exploration. As the current focus is mainly to clarify the role of chemical reactions on the mechanical responses rather than the other way around, the impact of mechanical rearrangement on chemical kinetics is not yet considered, and is beyond the purpose of this paper.

To capture the evolution of field variables (i.e. displacement, pore pressure and ionic concentration), the balance equations (equations (15) and (19)) need to be combined with the constitutive relations (i.e. stress–strain relationship in equation (25) and the adopted linear proportionality rules between Onsager's ‘forces’ and ‘velocities’ as stated in equations (46), (51) and (52)), and the loading and boundary conditions concerned. Consequently, the field equations could be formulated as a system of non-linear partial differential equations, which would require the development of a comprehensive numerical tool for solving coupled chemo-hydro-mechanical problems (Fernandez-Merodo et al., 2007; Tamagnini & Ciantia, 2016). However, this work mainly focuses on illustrating the model capacity to interpret chemo-mechanical coupled behaviour at the local level. Thus, to avoid numerical complication, the simulations in the following section fix the spatial fields by simplified homogeneous domains. Future numerical implementations are encouraged to study such processes by solving the full boundary value problem and explore the impact of spatial heterogeneity on the system response.

To demonstrate the capacity of the proposed modelling framework to capture the underlying physics of chemical dissolution, the proposed model is adopted to explain experimental observations using a series of debonding tests performed by Ciantia and his colleagues on calcarenite samples from the Mediterranean coast of Gravina di Puglia, Italy. In those tests the samples were exposed to a fully saturated acidic environment (Castellanza & Nova, 2004; Ciantia & Hueckel, 2013; Ciantia et al., 2015). The tested material is characterised by the presence of a network of calcite grains bonded by calcareous bridges, both of which contain a high percentage of calcium carbonate (CaCO₃) and are highly vulnerable to chemical dissolution when subjected to acidic flow. In the following, the proposed hydrodynamic model is adopted to capture the chemo-mechanical coupled behaviour of this calcarenite during chemical dissolution. Model parameters are calibrated based on uncoupled long-term debonding (LTD) tests, while coupled LTD tests are used to evaluate the model's performance. All the experimental data used in this section are taken from Ciantia et al. (2015) where detailed descriptions of the experimental set-up are presented.

The uncoupled LTD tests conducted by Ciantia et al. (2015) involved the immersion of a series of mechanically unloaded calcarenite samples in an acid reservoir with an initial pH value of 3·0. The samples were soaked for various times and the pH value of the reservoir was continuously monitored to follow the reaction rate. The weathered samples were then weighed after drying so that the mass of dissolved solid m_dis could be evaluated. The mass ratio of the dissolved solid and the initial reacting solid m⁰_s is used to reflect the degree of weathering as ξ_dis = m_dis/m_s⁰ (De Groot & Mazur, 2013; Ciantia et al., 2015). Accordingly, ξ_dis = 0 and ξ_dis = 1 correspond to intact and completely dissolved solid limits. Following this early stage, oedometric compression tests were performed on the dry, degraded samples to evaluate the residual strength at different degrees of weathering. This experimental procedure has adequately decoupled the mechanical and chemical processes, thus enabling to calibrate the corresponding parameters in isolation.

As the weathering stage of the above experiment is performed without mechanical loading, it can be used to calibrate the parameters controlling the chemical kinetics. Specifically, in the experiment additional acid was added in increments to the bulk fluid to maintain a strongly acidic regime (Ciantia et al., 2015). However, the frequency and exact amount per increment was not reported. Therefore, in the current simulation, instead of imposing additional acid surcharges, the initial volume of reservoir (V_R) is assumed to be big enough to provide sufficient hydrogen ions, while not letting external flux enter the system during weathering. Furthermore, for simplicity a spatial homogeneity is adopted in which the reservoir fluid is assumed to be always fully mixed without imposing any spatial gradients, thus avoiding having to solve the coupled diffusion–convection mechanism (i.e. J^β_i = 0 and J^γ_i = 0). As such, since the chemical reaction can only occur within the volume of the soil sample (V), and since the mass balance equation is used to capture the mass transfer within the reservoir, the mass balance in equation (19) should be corrected by the volume ratio between the soil sample and the reservoir (V/V_R)

As the molar quantity of the dissolved solid is proportional to the molar quantity of the consumed hydrogen ions by the stoichiometric number, the value of V/V_R can be back-calculated based on the mass ratio of the dissolved solid (ξ_dis) and the corresponding change in hydrogen ion concentration (Δx_H)

Experimental data from the weathering test (Ciantia et al., 2015) show that as the reservoir pH changed from 2·3 to 3·4 (i.e. Δx_H = 4·6 × 10⁻³ mol/l), 1·8% of the calcarenite sample (ϕ_s = 0·48) was dissolved (i.e. ξ_dis = 1·8%), so based on equation (55) V/V_R ≈ 1/200. The simulation in Fig. 4(a) suggests that b = 2 and k_c = 13 kg mol/(J cm³ h) are appropriate. Based on equations (52) and (53), the predicted dissolution rate (ξ) at pH = 2·5 (i.e. x_H = 3·2 × 10⁻³ mol/l) and room temperature is equal to 2·5 × 10⁻⁶ mol/(cm³ s), which is of the same magnitude as the calcite dissolution rate directly measured through experiments (Fredd & Fogler, 1998). The pH evolution in Fig. 4(a) implies that as the hydrogen ions are gradually consumed, the rising pH slows down the chemical dissolution, consequently bringing the pH value to a level (normally around 5·0) at which the chemical reaction becomes relatively slow due to the insufficient presence of the hydrogen ion. Within the time window of testing, it can even be viewed as the cessation of chemical dissolution. Furthermore, as suggested by the stoichiometric number of the chemical equation, every consumption of two H⁺ ions would generate one Ca²⁺ ion (Fig. 5(a)). Following equation (39), as no external flux enters the system during the soaking stage of the uncoupled LTD tests, and because the total number of ions lowers during the reaction, the osmotic pressure presents a decreasing pattern (Fig. 5(b)).

Based on the compression tests of the weathered samples, a Poisson's ratio of ν = 0·09 was deduced by Castellanza & Nova (2004), which is relatively low compared to the typical values reported for rocks, ranging from 0·1 to 0·4 (Gercek, 2007). Such underestimation might have arisen from the use of a soft oedometer, which enabled small radial deformation (rather than maintaining a strict K₀ condition) to allow measurement of the radial stress (Kolymbas & Bauer, 1993). For this reason, here the Poisson's ratio is set to be 0·25, as a representative value between the previously reported 0·1 and 0·4 values for rocks. The power coefficient n in equation (24) is fixed to 3 for calcarenite, as usually adopted for foams and soils (Gibson & Ashby, 1982; Rubin, 2001). Notice that a higher value of about n ≈ 10 would yield a much better fit, but this tuning is avoided, to ensure the model results in this paper would be viewed as predictions rather than as fits. However, the fact n ≈ 10 yields much better fits should encourage a careful experiment into the true value of n through elastic wave measurements on calcarenite samples. The other elastic properties are then determined based on the compression behaviour prior to elastic limits at various levels of weathering, as summarised in Table 2. Since the tested samples have a solid fraction of around 0·48, based on equation (26), K/K_s ≈ 0, which justifies the effectiveness of adopting the Terzaghi's effective stress principle to analyse the current experiments. Furthermore, according to the experimental data by Ciantia et al. (2015), the limit at which the intact sample under the oedometric conditions ceases to act as broadly elastic is reached when the axial stress σ_a approaches about 2700 kPa, which lowers to around 1800 kPa when ξ_dis = 6·2% (Ciantia et al., 2015). Since the proposed model behaves poroelastically, the model simulations in Fig. 4(b) focus on material behaviour prior to the elastic limit. These simulations show that, as weathering intensifies, the progressive dissolution of solid skeleton weakens the material strength and leads to the lowering of the compressive stiffness.

During the coupled oedometric LTD tests by Ciantia et al. (2015), the dry specimens were first loaded oedometrically till a prescribed axial stress (σ_a = 2000 kPa) where they were submerged by water flow to a fully saturated state under constant axial stress. Afterward, instead of water, a sustained acid flux was injected from the bottom to explore chemo-mechanical coupling effects. An input pH was imposed at 2·5, while the output pH was measured to be 5·0, at which due to the insufficient supply of hydrogen ion, the reaction rate was slow enough to be ignored within the experimental time frame. This infers that the pH spatially varies due to the presence of a transient acid flow, which would ideally require to treat the process as a boundary value problem instead of an elementary test (Ciantia et al., 2015). In other words, the process cannot be analysed solely at the constitutive level due to the lack of uniformity, but requires a proper numerical treatment to capture the spatial heterogeneity (Tamagnini & Ciantia, 2016). However, since the rate of the chemical reaction is much faster than the acid diffusion, the majority of the hydrogen ions are consumed by the material adjacent to the acid inflow (Fredd & Fogler, 1998; Tamagnini & Ciantia, 2016). It is thus anticipated that the acidity drops dramatically within a thin zone of intensive reaction while the pH values outside this reactive domain would approximately maintain constant as illustrated in Fig. 6 (Jyoti & Haese, 2021). Therefore, a model simplification can be made to entrap the whole process of chemical reaction within a localised basal layer susceptible to a constant pH, while letting the soil above remain as an intact solid block that does not experience chemical dissolution and follows the movement of the upper boundary of the localised layer due to compatibility (Fig. 6) (Chen & Buscarnera, 2022). However, owing to the lack of experimental evidence, the thickness of the reactive zone cannot be determined precisely. Instead, considering the spatiotemporal profile of pH simulated by Tamagnini & Ciantia (2016), the ratio of thickness between the sample height (H) and the localised zone (h) is approximated consistently by h/H ≈ 1/15.

In this context, flux through this reactive localised layer with volume V and surface S can be homogenised through the divergence theorem. Then, the mass balance equation in equation (19) can be adjusted for the finite reaction volume that is independent of gradient effects (⁠$∇iρ=∇iρβ=∇icγ=0$⁠)

where J_δ is the smeared flux of component δ entering and leaving the studied volume (positive for an outflow, negative for an inflow). Thus, combining equations (25), (33), (38), (56) and (69), the field equations for the finite volume could be obtained with further details shown in Appendix 4. With such simplifications, the coupled LTD experiment can still be simulated at the constitutive level while avoiding the computational complication associated with the spatial variation.

Besides, since the enforced axial stress (σ_a = 2000 kPa) is well below the elastic limit (≈ 2700 kPa at the intact condition, as discussed previously), it is reasonable to analyse the mechanical response as a poroelastic body. Admittedly, as the degradation exacerbates the material strength and diminishes the elastic limit, some inelastic processes might develop, such as material plasticity and grain crushing (Ciantia & Hueckel, 2013; Gajo et al., 2015). However, the impacts of these are neglected in the current work for simplicity, but future consideration is encouraged. Moreover, it is assumed that the material within the reactive localised zone experiences sufficient drainage without building up excessive pore pressures (i.e. P(t) = 0), and that the reaction-generated calcium ions can freely diffuse out of the system. However, no experimental measurement is available to assess the magnitude of the calcium ion flux, which requires further experimental exploration. Here, it is assumed that calcium ions can instantaneously escape from the system after generation (i.e. $c˙Ca=0$⁠). Otherwise, although not directly taking part in the chemical reaction, the ionic concentration of the anions (such as acetate derived from acetic acid) still needs to be accounted for, since besides hydrogen ions the acid injection brings in extra anions which later combine with the generated calcium ion to form an aqueous solution of calcium salt that would contribute to osmotic effects (Jyoti & Haese, 2021).

Considering the above-mentioned testing conditions, the simulation shown in Fig. 7 quantitatively predicts the measurements relatively well. As the acid progressively dissolves the solid skeleton, the sample deforms under constant axial load due to chemo-mechanical coupling and the radial stress grows to sustain the K₀ condition. Note that the simulation focuses on the axial strain developing only within the basal reaction zone of thickness h, while the upper part of the sample H − h is kept rigid. Therefore, the macroscopically measured axial strain (as reported in Fig. 7) is the one developing within the localised layer times the correction factor h/H.

In a similar way to the axial strain, the degree of dissolution ξ_dis also needs to be projected from the reaction zone to the overall sample, and Fig. 8(a) shows its evolution in which only the dissolved mass ratio at the end of the test is reported experimentally. Such predictions can certainly be improved by adjusting the reaction rate coefficient k_c to be state dependent, especially considering the impact of the surface area on the reaction kinetics, which, however, is left for future analyses. Fig. 8(b) implies that, as the acid is continuously depleted to dissolve the solid particles, the hydrogen ions need to flow into the system continuously to sustain a constant pH value (i.e. J_H < 0). Contrarily, since water and calcium ions are reaction products, they would freely escape from the system (i.e. J_w > 0 and $JCa>0$⁠) as governed by the prescribed boundary conditions, allowing instantaneous drainage and maintaining a constant calcium ion concentration.

This paper has introduced a thermodynamically consistent modelling framework to investigate the role of solid dissolution on the mechanical response of geomaterials by adopting a hydrodynamic procedure. Besides conventional state variables such as partial densities, momentum, elastic strain and entropy, the presence of chemical reaction required the incorporation of the mass fractions of ionic species in the pore fluid as additional state variables. The derivation showed that chemical reactions contribute to entropy production through ionic diffusion and chemical kinetics. To preserve the non-negativeness of the entropy production, simple yet effective linear proportionality relationships were adopted as the dissipative constitutive relations, which readily allows the recovery of some widely used empirical laws with clearly stated assumptions. It is shown that many of the said empirical laws tend to remain valid under the exact experimental conditions used for their original deductions. Beyond idealised conditions, it is instructive to favour physics-based models, thus highlighting the significance of the proposed model. For example, the widely recognised relationship in physico-chemistry to assess pressure dependency of chemical potentials is valid when the material is nearly incompressible, the commonly adopted Darcy's law is valid when the ionic species are nearly homogeneously distributed throughout the volume, and the standard Fick's law is valid when the pore fluid pressure does not show significant spatial gradients. Although the first scenario is typically true for geotechnical applications, pH and pore pressure profiles are not always homogeneous and thus require more generalised results. As such, the paper provides insights into cases where one normally does not (or technologically cannot reliably) obtain experimental data due to the unique boundary conditions.

Furthermore, the proposed theory holds generality, in that it should remain applicable as the basis for mechanically inelastic constitutive laws, such as elastoplasticity and breakage mechanics. Besides, the theory can be widely applied to investigate the chemo-mechanical coupling of any dissolution or precipitation process in porous systems, under both saturated and unsaturated conditions. However, the focus in this paper was on chemical degradation of saturated calcarenites induced by the dissolution of calcium carbonates in an acidic environment. Analyses were constrained to the elastic regime under a saturated condition, thus neglecting complications from mechanically inelastic phenomena such as solid rearrangement and grain breakage, as well as surface tension among solid–fluid–air interfaces. In so doing, the paper clarifies the significance of coupled chemo-mechanical processes and the role of chemical reaction on the mechanical response. The proposed model was used to interpret previously published experimental observations of LTD tests on calcarenite samples. Model parameters were calibrated based on the uncoupled LTD tests and the adopted model is shown to be able to capture the weakening of mechanical strength due to the loss of solid skeleton, and the decaying trend of chemical reaction due to the progressive consumption of acid. Furthermore, coupled LTD tests under drained conditions are simulated showing that the chemical reaction compresses the tested samples and tries to cause lateral expansion under unconfined conditions, thus leading to a rise of radial stress under an oedometric condition.

Although the consideration of only elasticity offers a satisfactory prediction of the experimental results, future extensions could look more deeply into the role of chemo-mechanical coupling on material inelasticity, including particle breakage and yielding response. Furthermore, instead of localising the response within a finite volume, numerical tools addressing the diffusive processes are desired to explore the spatiotemporal profiles of field variables, which would potentially benefit the analysis of various engineering applications and issues, such as landslides and sinkholes.

The authors would like to thank the Australian Research Council for funding through DP190103487, as well as Ebrahim Alaei, Leonardo Crespo and David Riley for fruitful discussions.

For systems undergoing multiple parallel chemical reactions, the mass conservation of each component involved is written as (De Groot & Mazur, 2013)

where M_δN_δjξ_j quantifies the density exchange rate of the δth component attributed to the jth chemical reaction; r represents the total number of parallel chemical reactions; M_δ is the molecular mass of the δth component; N_δj is the stoichiometric coefficient of the δth component in the jth chemical reaction, with negative values for reactants and positive values for products; and ξ_j (unit: mol/(m³ h)) indicates the rate of the jth chemical reaction, governed by the corresponding reaction kinetics, which is equal to 0 in the absence of chemical reaction. Typically, the dissipative density flux of the δth component is defined with respect to the barycentric motion

Combining the last two equations

Note that the mass of pore fluid is determined as the sum of all the ionic species involved and water, and therefore its change due to chemical reaction is (Bear, 1988)

Based on the definition of mass fraction of the ionic species in equation (6), it is found that ρ_γ = c_γρ_f. Combining this with equations (58) and (59), the evolution law of the mass fraction of the γth ionic species can be expressed as

Finally, the overall mass should be conserved and the total density follows

The first law of thermodynamics, which requires the overall energy to be conserved, is given by

where E_i is the energy flux. The gradient of energy flux can be derived by following the standard hydrodynamic procedure illustrated by Einav & Liu (2018) with additional state variables quantifying ionic diffusion and chemical reaction. Given the definition of the total strain rate $ε˙ij=−(∇ivj+∇jvi)/2$ and combining equations (8)–(21), (60) and (63)

where $A$ is the thermodynamic force

since the stoichiometric number N takes positive values for products and negative for reactants, by grouping the terms with the same sign of N, equation (65) essentially quantifies the sum of the total chemical potentials of products (μ_p) and reactants (μ_r), recalling that μ_p and μ_r take opposite signs, which is known as the driving force of the chemical reaction.

Then, the energy flux E_i is identified by grouping all the terms within the gradient, while all the remaining terms contribute to the rate of entropy production $R$

where the mechanical dissipation $D$ has been identified as

Based on the definitions in equations (9), (23), (27) and (28) the relationship between the intrinsic and the partial chemical potentials can be derived as

Furthermore, for saturated geomaterials the intrinsic pressure of all the involved phases should be identical based on the energy minimisation as proved by Jiang et al. (2017). Therefore, the thermodynamic pressure defined in equation (12) at the equilibrium state (i.e. v_i = 0) could be expressed as

According to equation (11) the total pressure p = p^e + 2u^e + P. As suggested by equations (24) and (25), the elastic energy density scales approximately as u^e ∝ (ε^e_ij)². Since for geomaterials, the elastic strains are normally small with the order O(ε^e_ij) ≪ 1, the contribution of the elastic energy density to the total pressure through thermodynamic pressure is relatively negligible compared to the contribution of the free energy through intrinsic pressure (Alaei et al., 2021). Thus, it is reasonable to replace u by u^m in both equations (9) and (12) so that

Combining the above with equations (4), (27) and (28) leads to

Based on equations (25) and (56) and considering only the elastic behaviour (i.e. ε^e_ij = ε_ij), the temporal evolution of the elastic stress is calculated as

Furthermore, equation (4) implies

Based on equations (6), (7) and (39) is alternatively expressed as

where $Kf=Kw+∑RTxγ$ is the bulk stiffness of the aqueous solution. For the dilute solution studied K_f ≈ K_w, while by inserting equations (36), (56), (71) and (74) into equation (73), the rate of change of the pore pressure becomes

where K_m = K_sK_f/(K_fϕ_s + K_sϕ_f) is the undrained bulk modulus of the solid–fluid mixture (Schmitt, 2015). After combining equations (56), (72) and (75), the field equations could be assembled in a matrix form

where the matrix coefficients are expressed as