Cracking and damage from crystallization in pores : Coupled chemo-hydro-mechanics and phase-field modeling

Cracking and damage from crystallization of minerals in pores center on a wide range of problems, from weathering and deterioration of structures to storage of CO2 via in situ carbonation. Here we develop a theoretical and computational framework for modeling these crystallization-induced deformation and fracture in fluid-infiltrated porous materials. Conservation laws are formulated for coupled chemo-hydro-mechanical processes in a multiphase material composed of the solid matrix, liquid solution, gas, and crystals. We then derive an expression for the effective stress tensor that is energy-conjugate to the strain rate of a porous material containing crystals growing in pores. This form of effective stress incorporates the excess pore pressure exerted by crystal growth – the crystallization pressure – which has been recognized as the direct cause of deformation and fracture during crystallization in pores. Continuum thermodynamics is further exploited to formalize a constitutive framework for porous media subject to crystal growth. The chemo-hydro-mechanical model is then coupled with a phase-field approach to fracture which enables simulation of complex fractures without explicitly tracking their geometry. For robust and efficient solution of the initial–boundary value problem at hand, we utilize a combination of finite element and finite volume methods and devise a block-partitioned preconditioning strategy. Through numerical examples we demonstrate the capability of the proposed modeling framework for simulating complex interactions among unsaturated flow, crystallization kinetics, and cracking in the solid matrix. c ⃝ 2018 Elsevier B.V. All rights reserved.


. Introduction
Growth of mineral crystals in pores can give rise to severe damage and cracks in the host material. ese coupled chemo-hydro-mechanical processes are now central to a number of problems in our society. A well-known example is weathering and deterioration of historic and building structures due to salt crystallization. Many of these structures are comprised of materials prone to invasion of salt water (e.g., stone), so they can be severely damaged when salt minerals grow inside the pores. See Fig. for example. Preventing this type of damage has been a critical element of conservation of cultural heritage and structures around the world [ -]. Also, crystallization of minerals in the subsurface can trigger ground heaving that severely damages buildings and geotechnical structures [ -]. Furthermore, reaction-driven cracking during mineral hydration, carbonation, and oxidation is a key consideration for deploying a promising strategy for geologic carbon storage that transforms CO into solid carbonate minerals [ ]. Addressing the problem of cracking and damage from crystallization in pores requires us, as a rst step, to understand its fundamental mechanism. Figure schematically shows how a crystal is present inside a pore space. Importantly, the crystal is usually con ned within a liquid solution, maintaining a thin liquid lm between its surface and the solid pore wall (see Scherer [ , ] for detailed explanations). e liquid solution surrounding the crystal allows it to continuously grow and push on the pore wall. is process of crystal growth generates an excess pressure on the solid matrix, which is commonly referred to as the crystallization pressure in a number of previous studies (e.g., [ -]). ese studies have proposed expressions of the crystallization pressure, showing that this pressure can far exceed the tensile strength of many porous materials. is a rms that the crystallization pressure is the direct cause of fracturing and damage during crystallization in pores, and thus, an accurate prediction of the crystallization pressure per se is an important research problem. However, an expression for the crystallization pressure alone is insu cient for tackling real-world problems such as those mentioned above, because the scales of these problems are orders of magnitude larger than the scale of pores. Also necessitated is a continuum-scale modeling framework that allows us to simulate and predict how crystallization pressures would evolve in space and time and ultimately a ect the problem at the eld scale.
e necessity of such a predictive modeling framework is the motivation of this work.
Continuum modeling of fracturing by in-pore crystallization poses two signi cant challenges. First, it requires a mathematical formulation that encapsulates complex interactions among chemical reactions, uid ow, and solid deformation in porous materials. A coupled formulation for such chemo-hydro-mechanical processes is not only di cult to develop in a theoretically consistent manner, but also challenging to solve numerically since the ow, transport, and reaction processes involve multiple length and time scales. Second, as can be seen from Fig. , the fracturing process of interest entails extremely complicated geometries that cannot be idealized as a set of sharp discontinuities. Obviously, algorithmic capture of such complex geometries are very unwieldy and onerous. Since these two types of challenges are interwoven herein, theoretical and computational modeling of crystallization-induced fracturing is a particularly demanding task.
Modeling frameworks that address both of these theoretical and computational challenges remain scarce. To our knowledge, Coussy [ ] introduced the rst chemo-hydro-mechanical framework for deformation and fracture from in-pore crystallization of minerals. His work presented signi cant contributions to theoretical aspects, but it remained unclear how the theory can be applied to construct an initial-boundary value problem of the crystallization-induced fracturing process. Some of more recent studies proposed computational models that paved the way to numerical simulation of the crystallization problem (e.g., [ , ]). Yet, challenges remain for both the theory and computation. First, theoretical formulations in past studies show signi cant disagreement, particularly with respect to the de nition of e ective stress which governs the constitutive behavior of the solid matrix. Second, the existing computational models appear insu cient for addressing complex cracking and damage processes from crystallization. For example, Koniorczyk and Gawin [ ] used linear elasticity without consideration of fracture. Derluyn et al. [ ] employed a simple local damage model along with elasticity, but their work was limited to identi cation of crack nucleation by comparing the e ective stress and the tensile strength in a -D setting. Such an approach would be inappropriate for delineation of the onset and evolution of damage and cracking zones in multi-dimensional problems. Also, their use of the standard, continuous nite element method (FEM) may not be an optimal choice, since it could su er from numerical stability problems arising from the advective transport of minerals. To address the stability problem, Derluyn et al. [ ] smoothed equations related to crystallization kinetics. However, they also found that the smoothing parameters can plague the physics signi cantly. Another issue is that their work used a one-way coupled staggered scheme, which is presumably because computational cost for fully coupled chemo-hydro-mechanics is prohibitively expensive without a carefully designed solution strategy. All of the aforementioned aspects indicate that a signi cant amount of more work is necessary to advance the theory and computation of crystallization-induced cracking and damage in porous materials. e purpose of this work is to develop a more theoretically grounded, computationally e cient modeling framework for crystallization-induced deformation and fracture in porous materials. To this end, we rst draw on recent advances in continuum modeling of coupled multiphysics in porous materials. For theoretically consistent modeling of the coupled chemo-hydro-mechanical problem at hand, we use continuum principles of thermodynamics which have been the basis of rigorous poromechanical frameworks in the literature [ -]. e use of thermodynamic arguments provides a systematic procedure to derive physically meaningful forms of e ective stress and constitutive laws for coupled multiphysics processes. Notably, the poromechanical frameworks developed in this way have demonstrated their ability to reproduce a variety of real-world observations (e.g., [ -]). en, for obtaining a stable numerical solution to the problem, we utilize the nite volume method (FVM) for the uid ow and transport problem which involves advection phenomena, and the standard nite element method for the solid deformation problem which does not. Furthermore, we introduce a three-eld block-partitioned solver that enables one to solve the coupled chemo-hydro-mechanical problem with an a ordable computational cost.
As for the modeling of cracking and damage, we adopt a phase-eld approach to fracture which has emerged as an e cient means for simulating complicated cracks without explicitly tracking their geometry (e.g., [ -]). is feature of phase-eld modeling is particularly desirable for our purpose since most (if not all) observed cracks due to in-pore crystallization are extremely di cult to delineate geometrically. Also, in most cases, crystallization in pores gives rise to distributed microcracks, resulting in rounded damage zones like those shown in Fig. . Such complex damage patterns can also be captured by the use of a phase-eld model, because it can be regarded as a particular class of nonlocal gradient damage models [ ]. As such, in this work we adopt a phase-eld model of fracture to simulate both distributed damages and localized cracks by in-pore crystallization, in a mesh insensitive manner. In doing so, we also propose a way to estimate the length regularization parameter of the phase-eld model, by drawing on an empirical relationship between the fracture toughness and the tensile strength of geomaterials. e paper is organized as follows. In Section , we formulate balance laws for a porous continuum containing liquid solution, gas, and crystals in the pores. Subsequently, in Section we use thermodynamic arguments in conjunction with the balance laws to derive suitable expressions for e ective stress and multiphysics constitutive relations. In Section , we derive a phase-eld formulation of brittle fracture as a balance law such that it can be augmented to the coupled chemohydro-mechanics formulation developed in previous sections. We then discretize the resulting formulation in Section via a combination of nite element and nite volume methods, and devise a block-partitioned iterative solver for the coupled problem. In Section we present numerical examples demonstrating the validity and performance of the proposed framework for modeling complex chemo-hydro-mechanical and fracturing processes from crystallization of minerals in pores.

. Conservation laws
In this section, we develop conservation laws for coupled chemo-hydro-mechanical processes in mineral-containing, partially saturated porous materials. We rst introduce a continuum mixture representation of this type of porous material, and then formulate balance laws for its mass, linear momentum, and energy. At this point, we note that the purpose of this and the next section is to formulate a general mathematical model amenable to being combined with various methods for modeling cracking and damage. Later in Section , we will adopt a speci c method, namely phaseeld modeling.

. . Continuum representation
Using mixture theory we conceptualize the material of interest as a multiphase continuum in which the solid matrix, liquid solution, gas, and mineral crystals are overlapped. See Fig. for an elementary volume representation of this four-phase mixture. We de ne the volume fractions of the constituent phases as where the index s refers to the solid matrix, l the liquid solution, g the gas, and c the crystals. It is noted that an index is used as a subscript when referring to an intrinsic property of a constituent phase, whereas it is used as a superscript when referring to a partial property of the mixture. We de ne the saturation ratios of the phases in the pore space-the liquid solution, the gas, and the crystals-as Let ρ s , ρ l , and ρ c denote the intrinsic mass densities of the solid, the liquid solution, and the mineral crystals, respectively. e partial mass densities of the constituent phases are given by where ρ is the mass density of the entire mixture.  Figure : Elementary volume representation of a four-phase mixture composed of the solid matrix, liquid solution, gas, and crystals. e liquid solution is also a sub-mixture of water and dissolved minerals.
As shown in Fig. , the liquid solution itself is also a mixture composed of water and dissolved minerals. Let the index w denote the water and d denote the dissolved minerals. eir global volume fractions are given by and their saturation ratios are We also de ne local volume fractions of the water and dissolved minerals as ( ) Similarly, their partial densities are de ned as Note that the intrinsic density of dissolved minerals ρ d is identical to ρ c , but we have used ρ d for notational purposes. Using the above variables, we can write the mass fraction of minerals dissolved in the solution as . . Balance of mass To derive conservation laws for this mixture, we use a kinematic description that traces the motion of the solid matrix. Let an overdot denote the material time derivative with respect to the motion of the solid matrix. en we can write balance of mass for the solid, the liquid solution, the gas, the dissolved minerals, and the mineral crystals aṡ respectively. Here, v is the velocity of the solid matrix, m α is the rate of mass production for phase α, and w α is the Eulerian relative mass ow vector of phase α with respect to the solid matrix, which is given by with v α being de ned as the velocity of phase α. For the dissolved minerals in the liquid solution, the relative mass ow is decomposed into convective and di usive/dispersive parts as where j is the di usive/dispersive mass ux of the minerals in the solution. At this point, we introduce a few assumptions that simplify mathematical expressions in the succeeding development. First, we assume that the water and crystal phases are incompressible and the solid matrix undergoes in nitesimal deformations. Second, adopting the assumption in the uid ow model of Castellazzi et al. [ ], we postulate that only the dissolved and crystallized minerals exchange masses among the ve constituent phases, i.e., m s = m g = and m l = m d = −m c . ird, we consider that momentum and pressure of the water and dissolved minerals in the liquid solution are under equilibrium.

. . Balance of linear momentum
Balance of linear momentum for the solid, liquid solution, gas, and crystals can be written as where σ α is the partial (Cauchy) stress tensor of phase α, a and a α is the acceleration of the solid matrix and phase α, respectively, and h α is the drag on phase α by the surrounding phases, which is subject to the constraint ( ) Summing up the above equations gives an expression for balance of linear momentum for the entire mixture, which reads Here, σ = σ s + σ l + σ g + σ c is the total stress tensor in the mixture.
. . Balance of energy Let P denote the total power, K the kinetic energy, and I the internal energy in the mixture. en balance of energy for the mixture can be expressed as ( ) e rate of change of kinetic energy in an arbitrary volume V of the mixture is given bẏ and the rate of change of internal energy isİ whereė is the rate of change of internal energy per unit total mass of the mixture. e total power is the sum of the mechanical power P m and the non-mechanical power P n , i.e., P = P m + P n . In this work, we consider isothermal conditions in which P n = . Let us introduce the in nitesimal strain tensor for the solid ε = ∇ s u = (∇ u + ∇ T u) where u is the displacement vector. Similarly we introduce in nitesimal strain tensors ε α for phase α. e mechanical power is then given by We now substitute above expressions intoİ = P −K, impose the balance of linear momentum for each phase, and localize the resulting integral. en, we can express the rate of change of internal energy as is energy balance equation will be the starting point of the next section where we derive suitable forms of the e ective stress tensor and constitutive laws.
. E ective stress and constitutive framework e purpose of this section is to identify physically meaningful constitutive relations that close the formulation. For this purpose, we begin by deriving a suitable form of the e ective stress tensor, which is a necessary step for constitutive modeling of a deformable porous material in ltrated by uid. We do this derivation based on a thermodynamic argument that the e ective stress should be energy-conjugate to the strain rate tensor. From the e ective stress derivation we identify several groups of energy-conjugate pairs, and we use them as a guide to develop constitutive relations for coupled chemo-hydro-mechanical modeling.

. . E ective stress
Our goal here is to nd the form of stress measure energy-conjugate to the strain rate tensor for the solid matrix, which, by de nition, corresponds to the e ective stress in a uid-in ltrated porous material. For this purpose we apply the procedure developed by Borja [ ], which was originally proposed for porous media without any crystals, to materials containing solid crystals in their pores. In doing so, we postulate that the phases in the pores-the liquid solution, gas, and crystals-can only carry volumetric stresses. is postulate is reasonable for most pore uids of interests. For the crystals, the postulate may be justi ed by the fact that they are usually oating in the liquid solution, as depicted in Fig. previously. See the same reasoning in Na and Sun [ ] for incorporating ice crystals into e ective stress in frozen soils. en, we can express the partial stress tensors of these in-pore phases as where p α is the intrinsic pressure of phase α and is the second-order identity tensor. Substituting the above expression into Eq. ( ) and neglecting thermal terms gives Eq. ( ), we can see that the second term on the right hand side is directly related to the solid deformation (because ∇ ⋅ṽ α = ∇ ⋅ v α − ∇ ⋅ v), whereas the last term is irrelevant to the solid deformation. us we now seek to extract the contribution of ϕ α (∇ ⋅ṽ α ) to the solid deformation. e mass balance equations for the liquid solution, gas, and crystals can be rewritten aṡ and rearranging the above equations gives ( ) Now we look for an alternative expression forφ α as it is also related to the solid deformation. Recalling that ϕ α = ( − ϕ s )S α for α = l , g, c, the material time derivative of ϕ α is given bẏ ( ) e remaining challenge is to have an explicit expression for the last term,φ s . Borja [ ] has shown that, in the case of barotropic (isothermal) ow for the solid, the mass balance equation for the solid can be expressed asφ where b = K K s , with K and K s being the bulk moduli of the solid matrix and the solid constituent, respectively. Inserting Eq. ( ) into Eq. ( ), we geṫ Using the above equations, we can rewrite Eq. ( ) as Here, the coe cients of ∇ ⋅ v can be simpli ed to where B = − b = − K K s is the well-known Biot coe cient. Finally, we can decompose ϕ α ∇ ⋅ṽ α in Eq. ( ) as Observe that the rst term on the right hand side contains ∇ ⋅ v, which is the rate of volume change of the solid matrix. Substituting Eq. ( ) into Eq. ( ), we can rewrite the energy balance equation as is the stress measure energy-conjugate to the strain rate tensor, which, by de nition, is the e ective stress. e "pore pressure" in this form of e ective stress, denoted byp, is which is the mean of the liquid, gas, and crystal pressures weighted by their volume fractions (saturation ratios).
In the absence of crystals (i.e., S c = ), Eq. ( ) boils down to the thermodynamically consistent e ective stress tensor in unsaturated porous media derived by Borja [ ], which specializes to the well-known Terzaghi's e ective stress [ ] in the limit of full saturation and B = . However, when a crystal is growing within a liquid solution in con ned pore space, the crystal growth exerts signi cant excess pressure on the solution lm between the crystal and surrounding solid particles.
is pressure-which is commonly referred to as the crystallization pressure in the literature-can make p c far greater than p l . en the e ective stress could become tensile even as the total stress is compressive, resulting in damage and fracture. Our derivation thus formally shows that the crystallization pressure can be a direct driver of deformation and fracture from crystallization in pores. Given its signi cance, the crystallization pressure is herea er denoted by p cr , i.e.,

. . Constitutive framework
Having derived the form of e ective stress, we further exploit the energy-conjugate pairs in the energy balance equation to identify suitable forms of constitutive laws. We particularly extend the procedure of Borja and co-workers [ , ], which has been used for unsaturated porous materials without crystals, to materials containing crystals in pores. is procedure simply interprets the type and form of constitutive relationship that each energy-conjugate pair suggests, without performing a thorough thermodynamic analysis. e reason is that the standard argument will eventually require variables in each energy-conjugate pair to be linked via a constitutive relation, for ensuring the second law of thermodynamics.
erefore, in what follows, we interpret the implication of each energy-conjugate pair, and adopt one of the widely used constitutive models in the literature. e selected constitutive model may not agree with the most general form allowed by a thermodynamic analysis. Note, however, that the speci c model is just chosen for constructing a particular class of the general framework we propose in this work-it can readily be replaced by another constitutive model depending on the purpose of modeling. Also, for brevity, the following discussion focuses on new aspects emerging from the presence of the crystals. Detailed explanations of the other aspects unrelated to the crystals can be found in [ -].
For a clearer interpretation, we rst simplify some terms in Eq. ( ). Since S l + S g + S c = , we haveṠ l = −Ṡ g −Ṡ c . us the second term in Eq. ( ) can be expressed as where the capillary pressure (suction) is de ned in the last term as p ca ∶= p g − p l . Also, the third and fourth terms in Eq. ( ) can be combined as erefore, under isothermal conditions, we can rewrite the energy balance equation as From the above equation we identify ve groups of energy-conjugate pairs, which means that ve types of constitutive relations are necessary to ensure non-negative entropy production. For each group, we can introduce a constitutive relation as follows. e rst energy-conjugate pair, from which we have de ned the e ective stress, is the mechanical power produced due to the deformation of the solid matrix. For this pair, it is natural to introduce a stress-strain relation of the formσ where C is a fourth-order tangent sti ness tensor. In this work, we shall assume that the solid behavior is linear elastic. is assumption may not allow us to fully describe the complex sti ness of geomaterials [ -], but it is o en good enough for modeling brittle fracture on which we focus in this work. e second pair contains the crystallization pressure and the rate of volume change of the crystals. Obviously this pair implies a constitutive law for the crystallization pressure. Notably, the derived energy-conjugate relationship is consistent with an equation for crystallization pressure put forward by Kelemen and Hirth [ ], which can be expressed as where ∆ψ is the change of Helmholz free energy and ∆V s is the di erence in volume between the solid products and the solid reactants. A number of speci c expressions have been proposed for the crystallization pressure (e.g., [ -]). Here we use the classical Correns' equation [ , ], given by where R is the gas constant, T is the temperature in Kelvin, V m is the molar volume of the mineral, and c eq is the equilibrium mass fraction at which the liquid is saturated by the mineral. is equation has later been shown to overlook the role of chemical activities in the crystallization process [ , ], but it gives reasonable predictions compared with experimental data of salt crystallization. e third pair contains the gas saturation ratio and the capillary pressure. Given that the gas saturation is determined by the liquid saturation, these variables may be related by a water retention law for unsaturated geomaterials. A common choice is the van Genuchten equation [ ], given by is equation requires four material parameters: the residual liquid saturation S l ≥ , the maximum liquid saturation S l ≤ , the scaling capillary pressure α ca , and the exponent n. Another exponent m is related to n via m = − n. is equation is originally developed for geomaterials in ltrated by freshwater, so it may be unable to accommodate some important aspects emerging from crystal growth in pores. However, the e ects of crystallization on the retention behavior have only become a subject of research recently (e.g., [ ]), and they are still far from being encapsulated into a water retention equation. For this reason, the original van Genuchten equation is used herein, but it can be modi ed easily when the e ects of in-pore crystals become quanti ed. On a related note, Derluyn e variables contained in the fourth group of energy-conjugate pairs are the relative velocity, pressure, and volume fractions of every phase inside the pores. is means that we need constitutive laws for ows of pore-lling phases relative to the solid matrix. For the liquid solution, we use the multiphase extension of Darcy's law, which can be written as where q l = ϕ lṽ l is the seepage velocity, µ l is the dynamic viscosity of the liquid solution, k r is the relative permeability, and k is the second-order absolute permeability tensor. We assume that the permeability tensor of the intact solid matrix is isotropic, i.e., k = k , and will discuss its anisotropic evolution by fracturing in the next section. Usually k is assumed to be constant when the solid deformation is in nitesimal. In our problem, however, crystals can clog much of the pore space.
To accommodate this clogging e ect, here we consider k a function of pore volume, adopting the Kozeny-Carman equation. e equation can be written in a normalized form where ϕ = − ϕ s is the porosity, and k and ϕ are the reference values of the absolute permeability and the porosity, respectively. e relative permeability is regarded as a function of saturation. When the water retention behavior is modeled by the van Genuchten equation, the relative permeability can be expressed as ( ) e relative ow of the dissolved mineral is comprised of convective and di usive/dispersive parts, as in Eq. ( ). e convective part is described by Darcy's law for the liquid solution presented above. For the di usive/dispersive part, we introduce a linear di usion equation In what follows, we shall assume that the gas pressure is atmospheric (i.e, p g ≈ ) and that the relative velocity of the crystal is negligibly small (i.e.,ṽ c ≈ ). ese assumptions, which have been introduced to other poromechanical models as well (e.g., [ , ]), allow us to neglect constitutive laws for the relative ows of the gas and crystal phases. Lastly, the h group of energy-conjugate pairs contains the rate of mass exchange, intrinsic density, pressure, and relative velocity of each phase inside the pores. Because m c = −m l and m g = , we only need to consider a constitutive law for m c . It is noted that the velocity has already been related to the pressure, and the pressure related to the saturation. is means that a constitutive law that relates the mass exchange, density, and saturation terms can satisfy this energy-conjugacy. In fact, by de nition, such a constitutive law corresponds to a kinetic equation for crystal growth and dissolution. Here we adopt the kinetic equation proposed by Espinosa-Marzal et al. [ ], which has later become the common choice of computational salt crystallization models (e.g., [ , , ]).
is equation can be expressed as where K c > and g c > are kinetic parameters, U ≥ is the supersaturation ratio, and U thr ≥ is the threshold value of U for crystal growth which depends on the type of mineral. Usually U thr ≥ for primary crystallization but U thr = once the crystallization process has begun. Multiple de nitions are possible for the supersaturation ratio U, and here we de ne U = c c to be consistent with Eq. ( ). Note that the rst equation in Eq. ( ) represents a crystal growth process (m c ≥ if U ≥ U thr ), whereas the second equation represents a crystal dissolution process (m c < if U < and S c > ). Note, however, that complex changes in the material's internal structure by dissolution (e.g., the formation of a sensitive clay by leaching) are beyond the modeling capacities of this formulation employing linear elasticity. So far, we have developed a general modeling framework for coupled chemo-hydro-mechanical processes in uid-in ltrated porous materials containing dissolved and crystallized minerals. In do-ing so, we have constructed a particular class of the framework by selecting a set of constitutive laws. It is again noted that our selection is just a speci c choice, and other sets of constitutive laws would work equally well. Similarly, while we adopt a phase-eld approach to fracture in the following, other approaches for similar purposes are also compatible with our development herein.

. Phase-eld formulation for fracture
In this section, we present a phase-eld model of fracture driven by the e ective stress derived in the previous section. To be consistent with the foregoing continuum mechanics approach, here we derive the phase-eld model as a balance law of microforces. Microforce balance derivations of phase-eld models of fractures have been presented in several previous studies [ , , , , ]. Among them, we adopt the derivation procedure of Choo and Sun [ ], which di ers from other microforce derivations because it views crack growth as a thermodynamically irreversible process. e motivation and implication of this derivation are explained in detail in Choo and Sun [ ].
Without loss of generality, in this section we consider a "dry" porous solid in which the liquid, gas, and crystal phases are absent. is simpli cation is just to prevent proliferation of numerous terms unrelated to the fracturing process. Again, an important premise in continuum poromechanics is that all mechanical processes-including the fracturing process which we will describe as the evolution of the phase-eld variable-are driven by the e ective stress. is means that terms other than the e ective stress and its related ones will not a ect the formulation that follow.
us, for brevity, the terms unrelated to the deformation and fracture of the solid matrix are omitted in this section.

. . Phase-eld approximation of fracture surfaces
We rst introduce a phase-eld approximation of fracture geometry, which, in essence, is an approximation of a sharp discontinuity as a di use interface. Let Γ denote a set of discontinuous fractures inside the body Ω. e total area of the fracture surfaces is given by which is an area integral over Γ. Calculating this integral during the evolution of fracture is an infeasible task because tracing the change of Γ is extremely di cult in most cases. To circumvent this di culty, we seek to transform an area integral over the evolving domain Γ into a volume integral over the xed domain Ω. For this purpose we de ne a phase-eld variable d ∈ [ , ], which denotes an intact state by d = and a fully cracked state by d = . Naturally d can be understood as a damage variable. Using this phase-eld variable, we introduce a crack density functional where the last integral is de ned over the volume Ω. In this work, we adopt a widely used crack density functional of the form where l > is the length parameter for the phase-eld regularization. e smaller the length parameter l, the closer the phase-eld approximation to the original sharp discontinuity.

. . Balance law derivation of phase-eld evolution
To derive a governing equation for the phase-eld variable, we make use of the microforce approach developed by Gurtin [ ], which has proven useful to derive phase-eld models [ , , , , ] and other types of models [ , ] within the framework of continuum mechanics. e rst step of this approach is to postulate the existence of a microforce system in which the phase-eld variable d (referred to as the order parameter in [ ]) is energy-conjugate to an internal microforce π and a surface microforce ζ. Consider an arbitrary volume V with boundary ∂V in this microforce system. e balance of microforce over the volume V is given by ( ) Denoting the unit normal vector of ∂V by n, we also introduce a microforce traction vector ξ such that ζ = ξ ⋅ n. en, by applying the divergence theorem and noting the arbitrariness of V, we obtain a localized form of Eq. ( ) as ( ) e internal and surface microforces, which are energy-conjugate to the phase-eld variable d, should evolve such that this balance law is satis ed. erefore this microforce balance law serves as a governing equation for the evolution of phase eld-equivalently, the damage and fracturing process. In the following, we derive speci c forms of π and ξ based on thermodynamic arguments. e mechanical power in the microforce system is given bỹ Incorporating this additional mechanical power as well as keeping the e ective stress power only, we can rewrite the balance of energy as We will exploit the second law of thermodynamics to derive expressions for the e ective stress tensor and the microforce variables. To this end, we rst consider a stored energy density function of the form Here, W(ε) denotes the strain energy stored in the undamaged material, and g(d) ∈ [ , ] is a so-called degradation function which should satisfy g( ) = and g( ) = . For now we consider general forms of W(ε) and g(d), and will discuss their speci c forms later in this section. At this point, it is noted that the energy used to create a fracture surface does not enter the stored energy function because crack growth is considered fully dissipative a priori herein. is is consistent with the stored energy functions used in variational frameworks for fracture (e.g., [ , ]), but di erent from those assumed in other balance law derivations of phase-eld models based on microforce arguments [ , , , ]. More speci cally, the stored energy functions in previous balance law derivations contain the fracture energy, so their derivations lead to a thermodynamic implication that crack growth is a reversible process in a rate-independent setting. is implication would be appropriate for clean fracture surfaces that can heal under highly controlled conditions. However, here we prefer to view crack growth as an irreversible process, since crack healing phenomena in geomaterials are usually inconsistent with the thermodynamic de nition of a reversible process. See Choo and Sun [ ] for a more detailed discussion on this aspect.
Having de ned the stored energy density, we can write the dissipation inequality as e time derivative of the stored energy function is given bẏ we get an alternative expression for the dissipation inequality as To ensure non-negative dissipation of the stress power irrespective ofε, the e ective stress should be related to the stored energy. is standard argument leads to a hyperelastic relation of the form Next, as done in [ , ] for deriving other types of models, we assume that the internal microforce is additively decomposed into two parts, the energetic (non-dissipative) part π en and the dissipative part π dis . In other words, π = π en + π dis . e same argument that we used to get Eq. ( ) yields the following expression for the energetic part:

Inserting Eqs. ( ) and ( ) into Eq. ( ) gives the reduced dissipation inequality of the form
Because the material is considered elastic, the dissipation is solely attributed to the evolution of the phase-eld variable d, or crack growth.
To arrive at speci c expressions for ξ and π dis , we now postulate that cracks are created in a way that they maximize the energy dissipation. is postulate is consistent with the Gri th theory of brittle fracture [ ], and has been central to several variational frameworks for phase-eld fracture (e.g., [ , ]). Our task is then to nd expressions for the microforce variables that maximize D f de ned in Eq. ( ). Equivalently, we seek to minimize the negative of the reduced dissipation functional, given by ( ) e arguments of this functional,ḋ and ∇ḋ, are indeed subject to a constraint in that these variables should form a phase-eld approximation of fracture. e speci c expression for this constraint can be obtained by taking time derivative of the crack density functional, Eq. ( ). is procedure giveṡ . Given this constraint, we can write a Lagrangian for this constrained minimization problem as where λ is the Lagrange multiplier of this problem. Invoking the stationary condition of this Lagrangian gives expressions for ξ and π dis as follows: ( ) e remaining task is to identify the meaning of the Lagrangian multiplier, λ, in the context of our problem. To this end, we substitute the results of Eqs. ( ) and ( ) into D f in Eq. ( ), which is the energy dissipation per unit volume. Integrating the resulting dissipation density over the domain Ω with the crack surface Γ d gives where the last approximation is attributed to the phase-eld regularization of discontinuous surfaces. Equation ( ) shows that the Lagrange multiplier λ can be interpreted as the energy dissipated by the creation of unit crack surface area. By de nition, this energy corresponds to the critical fracture energy in fracture mechanics. Let G c denote the critical fracture energy. Now we can express the internal microforce as and the microforce traction vector as Substituting these expressions into Eq. ( ), we can rewrite the governing equation for the phase- Notably, this equation is the same as the governing equation for a phase-eld model of brittle, rateindependent fracture obtained by variational and other balance law approaches, see [ , , , ] for example.
. . Stored energy function Now we consider speci c expressions for the stored energy function, ψ(ε, d) = g(d)W(ε). e stored energy function is central to the phase-eld model of fracture since it determines the e ective stress tensor σ ′ , as in Eq. ( ), and the energetic force driving the evolution of the phaseeld variable d, as in Eq. ( ). First, as for the degradation function g(d), we adopt the form most widely used by the phaseeld modeling community, given by ( ) It is noted that other forms of degradation functions have also been suggested, e.g., the cubic degradation proposed by [ ].
Next, we need to determine a suitable form of W(ε), the strain energy function decoupled from the phase-eld variable. In doing so, we should take into account that for physically realistic results, the strain energy associated with pure compression should not give rise to fracturing. For this reason, previous studies have proposed to decompose the stored energy function as W(ε) = W + (ε) + W − (ε), where W + (ε) is the fracturing part related to the tensile strain energy, and W − (ε) is the non-fracturing part related to the compressive strain energy. is decomposition has been mainly done via either of the following two schemes: one that uses the sign of principal strains Here we use the former one, and decompose the stored energy function of a linear elastic material as where ⟨⋅⟩ ± = (⋅ ± ⋅ ) ,λ andμ are the Lamé parameters, and ε a are the principal strains. A er decomposing the strain energy as above, we apply the degradation function to the fracturing part only, i.e., ( ) en the e ective stress is expressed as Likewise, the energetic microforce is expressed as Note that π en does not take the non-fracturing part of the strain energy, W − (ε).
. . ermodynamic restriction and crack irreversibility Our derivation leads to an expression for the dissipation inequality which needs to be satis ed for thermodynamic consistency. Rewriting Eq. ( ), we can express the rate of energy dissipation in a unit volume as which means that crack growth should be irreversible, as we postulated in the beginning of this derivation. Integrating this equation over the entire domain leads to an expression that is identical to the crack irreversibility condition presented in variational frameworks for phase-eld fracture (e.g., Eq. . us it can be concluded that our balance law derivation is consistent with the variational derivation with respect to the thermodynamic implication as well as the governing equation. Eq. ( ) can also be written as where we use an alternative form of the time derivative of the nonlocal term in the crack density functional. e above equation shows that the crack irreversibility condition boils down toḋ ≥ . To enforce this crack irreversibility condition,ḋ ≥ , we adopt the approach proposed by Miehe et al. [ ]. e approach is to make the energetic force driving the phase-eld evolution, which is denoted by π en in our derivation, a non-decreasing function even as W + (ε) is decreasing. Following this idea, we introduce a strain energy history functional H ≥ subject to the Karush-Kuhn-Tucker condition of ( ) Simply speaking, H is the maximum of the fracturing part of the strain energy during the course of loading. Replacing the stored energy term in Eq. ( ) with H gives a modi ed phase-eld equation of the form ( ) Eq. ( ) will be used as the governing equation for the phase-eld model in the sequel.
. . Permeability evolution by phase-eld fracture Phase-eld modeling of fracture in uid-in ltrated porous media should take into account the impact of the fracturing process on uid ow. Several types of approaches have been proposed for this purpose but with di erent emphases (e.g., [ , , , , , , ]). Most of them are concerned with hydraulic fracturing whereby the rate of mass transfer between the fracture and matrix systems (leak-o ) is o en far lower than the uid injection rate. e problem at hand, however, involves a quite di erent situation in that here crystal growth in matrix pores slowly drives fractures and signi cant uid ow between the fracture and matrix systems is expected. In addition, we face a new problem that growing crystals increasingly clog the fracture aperture. us, an approach that views the fracture and matrix pores as a whole-rather than distinguished domains described by separate governing equations-would be more appropriate for our purpose.
Given these considerations, here we take an approach that augments an anisotropic permeability tensor describing Poiseulle ow to the absolute permeability tensor in Eq. ( ). Such an approach has been advocated by Miehe and Mauthe [ , ] and Wang and Sun [ ], among others. We now write the absolute permeability tensor as Here, k matrix is the isotropic permeability tensor of the solid matrix which we have discussed in the previous section. On the other hand, k frac is anisotropic, and given by where n = ∇ d ∇ d is the unit normal vector perpendicular to the fracture direction, and w is the hydraulic aperture. Note that the hydraulic aperture herein should be a function of the crystal fraction as well as the fracture width. For simplicity, we assume that the aperture can be given by ( − S c )w, wherew is the hydraulic aperture in the absence of crystals. Forw, we adopt an equation proposed by Miehe and Mauthe [ ], specializing it to in nitesimal deformation conditions. e resulting form of the hydraulic aperture w is given by where l is the characteristic length of a line element perpendicular to the fracture. For simplicity, we assign l to be the mesh size h as in Miehe and Mauthe [ ].
. . Estimation of phase-eld modeling parameters for geomaterials Lastly, we describe how the phase-eld modeling parameters are estimated in this work. A standard phase-eld model of fracture requires two parameters: ( ) the critical fracture energy G c , and ( ) the length parameter l. While the two parameters have di erent origins-the former is from fracture mechanics theory and the latter is from a phase-eld approximation of fracture-they together determine the fracturing response of a material [ , , , ]. In particular, the peak stress under tension, which is usually reckoned as the tensile strength, is controlled by the fracture energy, length parameter, and elasticity moduli. It is thus of signi cant importance to assign a suitable combination of the phase-eld modeling parameters. Among them, however, the length parameter is uneasy to estimate because it is not a physically measurable quantity.
Here we devise a way to estimate the length parameter for geomaterials, drawing on an empirical relationship between the tensile strength and the fracture toughness. Let σ t denote the tensile strength of a geomaterial and K Ic denote its mode I fracture toughness. Experimental studies have found that, for many geomaterials under quasi-static conditions, the two properties are well correlated as a linear function [ -], given by Here, the unit of l is meter because β is given by a square root meter. is equation allows us to estimate l from experimental data of σ t and K Ic . Once we have estimated l in this way and known the tensile strength and elasticity parameters, we can nd the corresponding value of G c from Eq. ( ).
Interestingly, Eq. ( ) suggests that the length parameter is related to the relationship between the tensile strength and the mode I fracture toughness, rather than the values of these properties. Given that the linear relationship ( ) has been found for many types of geomaterials (though not necessary valid for all types of geomaterials, see [ ]), Eq. ( ) may serve as a useful guide to set a length parameter for phase-eld modeling of fracture in geomechanical problems.

. Discrete formulation
is section develops a discrete formulation for numerical solution of the chemo-hydro-mechanics and phase-eld model described so far. We begin by stating the governing equations and the relevant initial-boundary value problem. For robust and accurate spatial discretization of the coupled problem, we use the nite element method for the deformation and fracture problem (momentum balance and phase-eld equations), and the nite volume method for the ow and transport problem (mass balance equations). e resulting discrete system is solved by a staggered scheme that sequentially updates the chemo-hydro-mechanics problem and the phase-eld problem. To facilitate monolithic solution of the coupled chemo-hydro-mechanics problem, we design a blockpartitioned preconditioner which signi cantly speeds up iterative linear solvers.
. . Governing equations e coupled problem at hand is furnished by ve governing equations that follow. e rst one is the balance equation for linear momentum in the mixture. For simplicity we shall neglect the e ects of inertial forces and crystallization kinetics on the linear momentum. Substituting the expression for e ective stress in Eq. ( ) into Eq. ( ), we express the balance of linear momentum ( ) e second governing equation is the mass balance for the liquid solution, which was originally given by Eq. ( ). Since ρ l = ( − ϕ s )S l ρ l , we can expand this equation as Here, the second product can be simpli ed using the Biot coe cient B. To wit, we can rewrite Eq. ( ) asφ As for the primary variables of these governing equations, we select the following ve elds: the displacement vector of the solid matrix u, the pore pressure of the liquid solution p ∶= p l , the mass fraction of the dissolved minerals c, the saturation ratio of the crystals s ∶= S c , and the phase-eld variable d.
We complete the statement of the problem by prescribing boundary and initial conditions as follows. Let Ω denote the domain of interest and ∂Ω denote its boundary. e boundary is suitably decomposed as: displacement and traction boundaries, ∂Ω u and ∂Ω t ; pressure and liquid ux boundaries, ∂Ω p and ∂Ω q ; and mass fraction and mineral ux boundaries, ∂Ω c and ∂Ω j . ese decomposed boundaries satisfy ∂Ω = ∂Ω u ∪ ∂Ω t = ∂Ω p ∪ ∂Ω q = ∂Ω c ∪ ∂Ω j and ∅ = ∂Ω u ∩ ∂Ω t = ∂Ω p ∩ ∂Ω q = ∂Ω c ∩ ∂Ω j . e boundary conditions are given by where n is the unit outward normal vector and the hats denote the prescribed boundary values. e initial conditions are given by for all position vectors x ∈ Ω at time t = . For s and d, we consider no ux boundary conditions and zero initial conditions throughout.

. . Time discretization
We begin the discretization process by approximating the time derivatives in the residuals of mass balance for liquid, dissolved mineral, and crystals. Given that the chemical, hydrological, and mechanical processes in our problem may involve a variety of time scales, we use an unconditionally stable, rst-order backward Euler method. Consider a time increment ∆t from t n to t n+ . We discretize the time derivatives of the primary variables as v =u = u n+ − u n ∆t ,ṗ = p n+ − p n ∆t ,ċ = c n+ − c n ∆t ,Ṡ c = S c n+ − S c n ∆t , ( ) and the time derivatives of liquid saturation and liquid density asṠ l = (S l n+ − S l n ) ∆t andρ l = (ρ l ,n+ − ρ l ,n ) ∆t. All other variables are evaluated at time n + . For notational simplicity, herea er we drop the subscript n + for quantities at time t n+ .

. . Space discretization
For space discretization we use the nite element method for the deformation and fracture problem (momentum balance and phase-eld equations) and the nite volume method for the ow and transport problem (mass balance equations). is combination is motivated by mathematical natures of the two problems. e deformation and fracture problem is described by elliptic equations, so its discretization by the nite element method is an optimal choice. On the other hand, the ow and transport problem involves hyperbolic systems, for which the nite volume method is usually more appropriate than the (continuous) nite element method. is is mainly because nite volume discretization is robust in the presence of sharp gradients in the solution elds as well as locally mass conservative at the element level.
We use a single mesh for both nite element and nite volume discretization, as shown in Fig. . e degrees of freedom for the deformation and fracture unknowns-the displacement vector and the phase eld-are located at the element nodes, whereas those for the ow and transport unknowns-the pressure, mass fraction, and crystal saturation-are located at the element center.

Finite element degrees of freedom (displacement vector/phase field)
Finite volume degrees of freedom (pressure/mass fraction/crystal saturation) To begin nite element discretization of the deformation and fracture problem, we de ne the spaces of the trial solutions for u and d as where H is the Sobolev space of order one. Accordingly, the spaces of the weighting functions are de ned as rough the standard weighted residual procedure, we can readily develop the variational equations of the deformation and fracture problem as Here, we have presented the variational equations as residuals to solve these equations via Newton's method later. We then perform nite element discretization of Eqs. ( ) and ( ), and obtain the discrete residual vectors by assembling element contributions. e contributions of element e are given by where i denotes a shape function index. We use standard linear shape functions in this work. Next, we discretize the ow and transport problem via the nite volume method, on the same mesh used in the nite element discretization. We integrate the mass balance equations over each element, apply the divergence theorem to the uid ux terms, and multiply the residuals by the time increment ∆t. Element-wise contributions of the discrete mass residuals are then obtained as ∫ Ω e cS l ρ l B ∇ ⋅ (u − u n ) dV − ∆t ∫ Ω e m l dV + ∆t ∫ ∂Ω e (cw l ) ⋅ n e dA + ∆t ∫ ∂Ω e j ⋅ n e dA − ∆t ∫ ∂Ω e jĵ dA = , ( ) where n e is the outward normal to the boundary of element e. As depicted in Fig. , the interpolation functions for the pressure, mass fraction, and crystal saturation elds take a constant value of over element e and at all other elements. us the volume integrals can be evaluated much like integrating nite elements with a piecewise constant shape function. e surface ux integrals are evaluated as a sum of interelement uxes between element e and its neighboring elements f . Let n face denote the number of faces of element e and n e f denote the outward normal vector at the interface ∂Ω e f . e integral of the liquid mass ux is then expressed as ∫ ∂Ω e w l ⋅ n e dA = We apply a two-point ux approximation scheme for multiphase ow in porous media as described in [ , ]. Since w l is governed by Darcy's law, (w l ) e f is given by Here, Φ e = p e + (ρ l ) e gz e and Φ f = p f + (ρ l ) f gz f are the ow potentials at elements e and f , respectively, where g denotes the gravitational acceleration and z denotes the elevations. Γ e f is the transmissibility at the interface ∂Ω e f , which is multiplicatively decomposed into the geometric transmissibility T e f and the uid mobility ( where A e f is the area of the interface ∂Ω e f , l e and l f are the distances between the center of the interface and the center of elements e and f , respectively, and k e and k f are permeabilities of elements e and f , respectively. Note that this term is independent of the ow direction. On the other hand, the uid mobility (λ l ) e f is evaluated in an upstream weighting (upwinding) manner as follows: ( ) e remaining ux terms can be computed in a similar manner. For evaluating the surface integral of cw l ⋅ n e , c e or c f is multiplied to (λ l ) e f in Eq. ( ) depending on the ow direction. Modi cation of the above equations to calculate the interelement uxes of j = −ρ l D ∇ c is straightforward, and omitted for brevity.

. . Fully discrete form and staggered scheme for phase eld
We now develop the matrix form of the discrete residuals as a linear system that needs to be solved in each Newton iteration. In doing so, we split the overall system into two parts-the phaseeld part and the rest-and apply a staggered solution method developed for phase-eld modeling of fracture. e primary motivation of using a staggered scheme is its computational robustness: it is far more convergent than a monolithic scheme during the evolution of phase-eld fracture. Such robustness is particularly more desirable for phase-eld modeling of fracture in strongly coupled multiphysics problems [ , ]. e phase-eld system, which is uncoupled from other variables, is given by where J dd is the Jacobian matrix of R frac with respect to the phase-eld variable d, and ∆ denotes a Newton increment. Indeed the phase-eld system is linear, so d can be updated by solving Eq. ( ) once. Next, the rest part-the chemo-hydro-mechanics system-is given by the following block-partitioned system: where J (⋅)(⋅) denotes a sub-matrix in the Jacobian matrix. Speci c expressions for the sub-matrices are straightforward to obtain and omitted for brevity. is chemo-hydro-mechanics system is usually nonlinear due to water retention characteristics and crystallization kinetics, among other reasons. It is noted that J s(⋅) and J (⋅)s are nonzero only when the crystals are growing or dissolving in the pores. Extending the staggered scheme proposed by Miehe et al. [ ], we proceed the numerical solution from time t n to t n+ through the following three sub-steps: . Determine H in the phase-eld system with the displacement variable u at t n .
. With this H, update the phase-eld variable d at t n+ , by solving Eq. ( ).
. With the updated d, update all other variables u p c s at time t n+ , by Newton's method solving Eq. ( ) at each iteration.
In the last step, we solve the coupled chemo-hydro-mechanics system in a monolithic manner that updates the primary variables in Eq. ( ) together in each iteration. However, other methods that sequentially solve the coupled system (e.g., [ ]) can also be used for the same purpose. Such sequentially-implicit methods may be preferred particularly when separate code is available for each physics, e.g., one for the deformation and fracture problem and another for the ow and transport problem. Yet, the analysis of White et al. [ ] suggests that, with respect to computation cost per se, a sequential solution approach is suboptimal to a monolithic one employing a decent block-partitioned preconditioner. is motivates the following discussion on a block-partitioned preconditioner for our problem at hand.

. . Block-partitioned preconditioner
We now seek to design a quality preconditioner for the Jacobian matrix in the coupled chemohydro-mechanics system given by Eq. ( ). To this end, we extend a block-partitioning strategy originally developed for poromechanics [ -] to our problem which involves more coupled physics. For the purpose of preconditioning, we rst consider a subset of the problem whereby crystal kinetics are insigni cant (i.e., J s(⋅) → and J (⋅)s → ). In such cases, the Jacobian matrix reduces to which is a by block-partitioned matrix. Our particular interest is in developing a block lowertriangular preconditioner of the form which, when multiplied to the Jacobian matrix J, results in a matrix very close to the block uppertriangular matrix obtained from the LDU factorization of J. In other words, we want a block lowertriangular matrix P − such that irrespective of U , U , and U . If this equation is satis ed exactly, the preconditioned Jacobian has a single eigenvalue of and the Jacobian system can be solved by a couple of Krylov iterations. Such an exact preconditioner, however, is generally very expensive to construct. erefore, to nd an e ective yet inexpensive preconditioner, we take the following three-step approach: . Find the block preconditioner that exactly satis es Eq. ( ).
. Approximate dense terms in the exact block preconditioner.
. Replace inverse operations in the approximated block preconditioner with their own preconditioners.
is approach has been put forward by White and co-workers [ -] for coupled poromechanics, and it has proven very useful for arriving at a scalable preconditioner. In the following we apply this approach to the coupled chemo-hydro-mechanics problem of our interest.
First, we solve for Eq. ( ) to nd a block lower-triangular matrix P − that exactly satis es this equation. We then obtain six equations for the six sub-matrices of P − as follows: ( ) e rst equation obviously leads to X = J − uu . From the second and third equations we get Here, J pp − J pu J − uu J up is the Schur complement of the block J uu of the upper le by block matrix, which will be denoted by S u in the following. Solving the rest, we get Similar to before, Z appears in L and L . Second, we approximate dense terms in the exact preconditioner derived above. Of particular interest is in nding good approximations of those containing J − uu multiplied by coupling matrices, e.g., S U = J pp −J pu J − uu J up . ese terms and their inverses appear all blocks of the exact P − except the X block, but their exact values are overly expensive to be used for a preconditioning purpose. us we will approximate these terms, drawing on the " xed-stress" split scheme for poromechanics which has been widely used as a sequential solution method [ , ] and recently rephrased as a blockpreconditioning approach [ ]. e split scheme goes as follows. If the mean volumetric stress is assumed to be xed during uid ow, the divergence of the solid velocity can be related to the change of pore pressure only, i.e.,σ Observe that the displacement vector u is now eliminated from the original mass balance equations. Discretization of these equations gives sparse approximations of the coupling terms containing J − uu . For example, the approximated Schur complement, denoted byS U , is given by where V e is the volume of elements assembled in an element-wise manner (corresponds to the mass matrix of the pressure shape functions in nite elements). Here, the original dense term J pu J − uu J up is approximated by a sparse term (the second term).
is approximation reduces to that derived by White et al. [ ] in case the material is fully saturated (S l = ) and ρ l is constant. Similarly, we approximate other terms as Accordingly, the Z block is also approximated as is term may be further simpli ed toJ − cc or (J cc −J cp diag(S − u )J pc ) − . We use the former one in this work.
ird, we replace inverse operations in the approximated block preconditioner-which are usually expensive to compute-with their own preconditioners. Speci cally, where P − (⋅) is a preconditioner for (⋅) − . We use algebraic multigrid preconditioners in this work, but a number of other combinations of sub-preconditioners can be e cient as well. Discussions on the choice of sub-preconditioners are presented in [ -] in the context of poromechanics problems.
Following the three-step procedure described above, we nally arrive at a speci c expression for Eq. ( ) as ( ) It would be illuminating to demonstrate how this block-preconditioner operates in real problems. Indeed, what is needed by an iterative solution method is the matrix-vector product of the preconditioner and an input vector x, say P − x = y. is operation can be written as ( ) e vector y can be e ciently updated via the following three steps: .
Notice that, from the second step, we use the block(s) of y that we updated in previous step(s)-y u in the second step, and y u and y p in the third step. In this sense, this block-preconditioning approach can also be interpreted as a sequential solution method. See White et al. [ ] for a thorough discussion on this interpretation. e repeated pattern in this sequence can be further applied to update the additional block related to the crystal kinetics, which we have observed satisfactory results in the numerical examples that follow. e proposed block-preconditioner has shown decent performance (usually less than Krylov iterations for a tolerance of − ) for a number of chemo-hydro-mechanical problems we have tested, and it is believed to be a good choice for other problems sharing a similar mathematical structure. However, much more extensive studies are needed to verify and improve the performance of the current preconditioner-especially with respect to its scalability in large computing platforms, as done by [ -] for poromechanical problems. To better focus on the scope of the current work, we leave this topic as a future research direction.

. Numerical examples
In this section, two numerical examples are presented to validate and demonstrate the capability of the developed computational model for simulating complex interactions among unsaturated ow, solid deformation, fracturing, and crystallization in pores. e rst example serves as a benchmark problem for verifying the numerical implementation and validating the mathematical model. We simulate a validation problem in unsaturated poromechanics, because a benchmark setting for chemo-hydro-mechanics involving in-pore crystallization is yet to emerge to the best of our knowledge. Subsequently, in the second example we incorporate the ow, transport, and crystallization of minerals and the resulting fracturing process in the host material. We particularly consider the problem of capillary in ltration of saturated solution into a porous medium, which has been commonly used by a number of experimental and numerical studies on crystallization of salts in geomaterials (e.g., [ , , ]). e main purpose of this example is to showcase the performance of the computational model for simulating the development of fully cracked regions ensuing crystallization in pores.
e numerical examples have been prepared using Geocentric, a massively parallel nite element code for geomechanics that has been used in a number of previous studies (e.g., [ , , , , , , , ] . . Drainage of freshwater in a porous column Our rst example has two purposes: ( ) to verify the numerical implementation of combinednite element and nite volume discretization in the current code, and ( ) to validate the mathematical model's capability for capturing coupled hydro-mechanical responses that can lead to crystallization of materials. For the latter purpose, we particularly focus on a drainage (saturation decreasing) process whereby supersaturation and crystallization can take place.
To achieve these two purposes simultaneously, we consider a benchmark problem of unsaturated poromechanics that simulates Liakopoulos' drainage experiment of a sand column [ ]. is problem has been used by a number of previous studies to validate their own combinations of governing equations, e ective stress, water retention, relative permeability, and others (e.g., [ -]). Here we use this problem in the same context. To verify the implementation, we simulate the drainage experiment with the current code using the nite element method for the solid deformation problem and the nite volume method for the uid ow problem, and compare the results with those obtained by unsaturated poromechanics code using mixed nite elements for both the solid deformation and uid ow problems (see [ , , , ] for details of the mixed nite element discretization). It is noted that, for problems like this one that do not involve a strong advection phenomenon, nite volume and nite element results should converge to the same solution. e numerical results are also compared with the experimental data of Liakopoulos [ ] to validate the mathematical model. Figure shows a schematic illustration of the problem setup. A . m wide and . m tall column is prepared such that it is initially under equilibrium with zero pore pressure. e test begins by making the top boundary undrained, which leads to drainage of the pore uid through the bottom boundary. During the test the bottom boundary is kept zero pressure and xed, whereas no ux and no lateral displacement are allowed along the side boundaries. Both of the uid ow and solid deformation in this problem are one-dimensional, and the column is discretized by quadrilateral elements along the height. e same mesh is used for both nite element/ nite volume and mixed nite elements simulations. With a uniform time increment ∆t = minute, we simulate the experiment until the drainage time reaches minutes. Material parameters are assigned close to the properties of the Del Monte sand tested by Liakopoulos [ ]. e material is composed of sand grains and fresh pore water whose intrinsic densities are ρ s = kg/m and ρ w = kg/m , respectively, with porosity ϕ = − ϕ s = .
. As for hydraulic parameters, we use the absolute permeability k = . × − m (isotropic), the dynamic viscosity µ l = − kPa⋅s, and the van Genuchten model parameters S = , S = , α ca = kPa, and n = . Because no mechanical data is available for the tested sand, we assume Young's modulus E = kPa and Poisson's ratio ν = . . Initially, the domain is fully saturated with zero pore pressure (p = ). Once the simulation begins, the pore uid is drained through the bottom outlet lter by gravitational force.
In Fig. we present simulation results obtained by the nite element/ nite volume and mixed nite element schemes, along with the experimental data of Liakopoulos [ ], in terms of time evolutions of pore pressure variation along the height and uid velocity at the bottom outlet lter. First, we nd that the two numerical schemes yield virtually identical results, thus verify our implementation of combined nite element/ nite volume discretization. Second, we see that the mathematical model well reproduces the spatio-temporal evolution of pore pressure measured in the experiment. In Fig. a the agreement between the simulation and experiments is rather qualitative in the beginning, but it becomes increasingly quantitative as time proceeds. Figure b shows that the calculated outlet velocities are in an excellent agreement with the measured data throughout the test. ese results demonstrate that, even when the solid behavior is grossly simpli ed, the computational model can capture salient physics of coupled solid deformation and unsaturated ow in real geomaterials.

. . Capillary in ltration of salt water and crystallization-induced cracking
We now proceed to the simulation of chemo-hydro-mechanical processes involving crystallization in pores and resulting damage in the solid matrix. Speci cally, we simulate a laboratory-scale problem analogous to a typical experiment studying growth of salt crystals in geomaterials [ , ], in which a specimen is in ltrated by salt water by capillarity and then damaged by crystallization of salt minerals in pores.
Figure illustrates the setup of this problem. We consider a porous rock column of width . m and height . m. e column is initially lled with freshwater and air, manifesting a uniform capillary pressure of kPa. To emulate capillary rise of a mineral solution, we prescribe the pressure and concentration at the bottom boundary as kPa and c eq (equilibrium solute mass fraction), respectively. e liquid solution is not allowed to ow outside through the lower half of the lateral boundaries as well as the top boundary, like in the experiment of Noiriel et al. [ ]. In contrast, the solution is subject to a constant outward ux boundary condition ofq = − × − kg/m /s in the upper half of the lateral boundaries. roughout the domain the dissolved minerals are subject to no ux boundary conditions such that the minerals stay in the pores during an out ow of the liquid solution.
ese boundary conditions are intended to drive supersaturation of the solution-and ultimately crystallization in pores-in the upper half of the specimen. Exploiting the symmetry of the problem, we model the le half of the domain and re ect the result in the post-processing stage. e half domain is discretized by , quadrilateral elements of uniform size (mesh diameter h = . cm), which leads to , degrees of freedom for the ve primary variables. We assign the material's properties as follows: solid density ρ s = kg/m , initial porosity is ϕ = − ϕ s = . , absolute permeability k = . × − m , and van Genuchten parameters S = . , S = . , α ca = kPa, and n = . . e linear elasticity parameters of the solid matrix are assumed to be E = GPa and ν = . . e properties of the pore water remain unchanged from the previous example. As for the mineral, we consider salt (sodium chloride), and adopt its parameters mainly from Castellazzi et al. [ ]. ey are: the mineral density, ρ m = kg/m , the molar volume, V m = cm /mol, the equilibrium solute mass fraction, c eq = .
[ ] and σ t = MPa, we get the phase-eld model parameters as G c = . J/m and l = . × − m. It is noted that the ratio of the length parameter to the mesh diameter l h is greater than , which has been shown to give su cient accuracy [ ].
Once the simulation begins, the salt water in ltrates the specimen from the bottom boundary. Figure shows the evolution of the supersaturation ratio U until minutes. At minutes and minutes the maximum supersaturation ratio equals unity, since no out ow is allowed through the lower half boundaries for both the solution and dissolved minerals. However, when the solution reaches the upper half boundaries, the water in the solution ows outside through the boundaries, whereas the dissolved minerals remain inside the specimen. is type of out ow makes the supersaturation ratio higher than unity. As a result, the contour at minutes show supersaturated zones in which U > . We can also nd that the supersaturation ratio exceeds U thr = . in some regions, which means that some minerals have been crystallized. Figures and present the crystal saturation ratio, S c , and the crystallization pressure, p cr , at three time instances from minutes to minutes. Figure shows that during this time crystals grow continuously at the upper half of the lateral boundaries, which may be viewed as salt e orescence. We can also see from Fig. that this crystal growth leads to an increase in crystallization pressure. Recall that the crystallization pressure enters the mean pore pressure in e ective stress. As such, this increase in the crystallization pressure can give rise to tensile e ective stress that drives the evolution of the phase-eld variable.
In Fig. we plot the evolution of the phase-eld variable d from minutes, with a uniform time interval of . minute. Until . minutes the phase-eld value does not reach . , so the material may be regarded as partially damaged by microcracks. A er . minute, the phase-eld value evolves to at some locations, which means that the material is fully cracked therein. e cracked zones enlarge further in the next . minute. It is noted that the shape of the cracked regions at . minutes resembles characteristic cave-like damage zones in salt weathering. To the best of our knowledge, this is the rst numerical simulation of the evolution of damage and fully-cracked zones by in-pore crystallization of minerals. e simulation is terminated at this point, because many assumptions in the constitutive models break down upon fracturing of the material. For example, the crystallization pressure equations in the literature have been proposed for a crystal con ned in small pore space such that a solution lm exists between the crystal and the solid matrix. erefore, when cracks develop signi cantly as . minutes in Fig. , these equations might become invalid or the crystallization pressure might disappear. Also, we have assumed that the kinetic parameter K c is constant, but it is expected to evolve by the fracturing process. Note, however, that all these limitations stem from the lack of constitutive relations for crystals within fractures, not from the modeling framework. Once constitutive models that address the aforementioned issue become available, they can be cast into the current modeling framework. Lastly, we would like to discuss computational aspects of the problem. As shown in Fig. , the phase-eld variable has evolved from partial damage to full cracks dramatically within the last . minute. To capture such an accelerated cracking process, we have used an adaptive time stepping algorithm proposed by Borden et al. [ ]. As a result, during the cracking process the time increments have been reduced to very small numbers (orders of magnitude smaller than a minute), giving rise to undrained deformations in which the relative ow of the pore uid is negligible. is aspect further justi es our choice of nite volume discretization for the ow problem, because mixednite elements for coupled poromechanics are subject to an inf-sup stability condition in undrained conditions [ , -]. We also would like to mention that solving a large number of steps in the crack development stage has been made a ordable thanks to the block-partitioned preconditioner described in the previous section. is preconditioner has allowed us to solve each linear system via a couple of dozens Krylov iterations until large cracked zones emerge.
. Closure is paper has presented a theoretical and computational framework for modeling cracking and damage in porous materials by in-pore crystallization of minerals. e framework combines a chemo-hydro-mechanics approach with a phase-eld description of fracture. Particular contributions of this work include: ( ) derivation of a thermodynamically consistent e ective stress tensor in porous materials containing growing crystals, ( ) identi cation of state variables that must be linked via constitutive laws, and ( ) a block-preconditioned iterative solver that facilitates numerical solution of the fully coupled chemo-hydro-mechanics equations. e computational model has been demonstrated to be capable of simulating the onset and evolution of cracking and damage from in-pore crystallization of minerals.
In this work we have proposed a general modeling framework and then constructed its particular class by employing relatively simple constitutive models of crystallization kinetics and solid deformation. We would like to note that use of more advanced constitutive models for these physics may lead to other classes that have better predictive capabilities. Particularly, if quanti able, porescale characteristics need to be incorporated into the crystallization kinetics model. An important example is the pore size distribution because crystallization in real geomaterials begins from smaller pores. e orientation of the mineral crystal may also exert control on this problem [ ]. As for the solid model, incorporating pressure-dependent plasticity will allow us to accommodate the e ect of con ning pressure on fracturing [ ], which is important for deep subsurface problems. Beyond constitutive models, the temperature eld must be augmented to the formulation when heat ow is signi cant. For problems involving complex drying phenomena, one needs to introduce additional phases for vapors and models of moisture transport [ -]. e framework developed in this work is amenable to incorporating these more complicated physics. ---. Support from these sponsors is greatly appreciated. e views and conclusions contained in this document are those of the authors, and should not be interpreted as representing the o cial policies, either expressed or implied, of the sponsors, including the Army Research Laboratory or the U.S. Government. e U.S. Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation herein.  [ ] S. J. Kowalski, ermomechanics of drying processes, Springer, .