Joumana Ghoussoub

Yves M Leroy
email: leroyy@lms.polytechnique.fr

Solid-uid phase transformation within grain boundaries during compaction by pressure solution

Keywords: A. Grain boundaries, A. Phase transformation, B. Geological material, B. Finite elements, C. Stability

The overall compaction of porous rocks due to intergranular pressure solution (IPS) results from the dissolution of minerals within contact regions and the di usive transport through the grain boundary of the dissolved species towards the uid-ÿlled pore space. The grain boundary structure can be imagined to be composed of dry contact zones, thin uid ÿlms and uid-ÿlled cavities. The connectiveness and tortuosity of this structure determine the e ective di usivity of grain contacts and thus the potential of porous rock to compact by the action of IPS. The evolution in time of the grain-boundary structure, and thus of the e ective di usivity, is discussed here with the help of two 2D initial-and boundary-value problems which are solved by analytical and numerical means. The evolution of the solid-uid interfaces within the grain boundary is governed by a phase transformation between the non-hydrostatically stressed elastic solid and the trapped uid assumed in mechanical equilibrium. The characteristic time is provided by a linear kinetic law. The evolution of the structure away from a state of thermodynamic equilibrium during a loading normal to the grain boundary is found to occur in two steps. The ÿrst one consists of a di use morphology evolution in time and results in an enhancement of any initial stress concentration. The second step is characterized by a rapid and localized dissolution in the region of stress concentration. The latency period prior to localization is governed by the magnitude of the non-hydrostatic remote stress as well as the microstructural geometric factor responsible for the initial stress concentration at the solid-uid interface. The localized dissolution is shown to provide a mechanism for the uid to penetrate a previously dry contact region by marginal dissolution and thus to create a uid ÿlm. However, the newly formed thin uid layer is found to be unstable pointing to a possible repeated reorganization or dynamic evolution of the grain boundary internal structure during the action of IPS.

Introduction

The objective of this paper is to study the evolution with time of the grain boundary internal structure under conditions typical of the action of intergranular pressure solution (IPS).

The deformation mechanism IPS di ers from classical Coble creep mainly by the presence of a uid phase within the grain boundary. IPS consists of the dissolution of minerals at stressed contact points, transport of the dissolved material through the grain boundary and its redeposition elsewhere. It occurs in the Earth's crust at geological rates and results in lithiÿcation and compaction of low grade metamorphic rocks [START_REF] Rutter | Pressure solution in nature, theory and experiment[END_REF]. IPS is partly responsible for the reduction in porosity of sedimentary rocks, a key parameter in determining the quality of potential hydrocarbon reservoirs. IPS has also been invoked to explain the creeping and the strength recovery of faults [START_REF] Sleep | Ductile creep, compaction and rate and state dependent friction within major faults[END_REF]. Finally, it is active in halite [START_REF] Urai | Weakening of rocksalt by water during long-term creep[END_REF], a material which has been studied extensively for the safety of nuclear and industrial waste disposals.

Rheological models for IPS often include the assumption that the grain boundary is a Gibbs interface having some e ective di usivity responsible for the transport of dissolved minerals [START_REF] Rutter | The kinetics of rock deformation by pressure solution[END_REF]Lehner and Bataille, 1984/85;Spiers and Schutjens, 1990;[START_REF] Paterson | A theory for granular ow accommodated by material transfer via intergranular uid[END_REF]. The di usivity coe cient is a phenomenological constant and is deduced from experiments performed on a collection of compacting grains. This choice of a constant value cannot model the exact nature of the action of IPS at the grain boundary scale and the potential healing of grain boundaries which would prevent the action of IPS. Three models have been proposed to describe the action of IPS at the grain boundary scale. The ÿrst one assumes that the dissolution proceeds from the pore space towards the inner contact region [START_REF] Bathurst | Diagenetic fabrics in some British Dimantian limestones[END_REF] and could be coupled with some crystal plasticity. The key feature of the second model is the existence of an adsorbed uid ÿlm at the grain boundary [START_REF] Weyl | Pressure solution and force for crystallization-a phenomenological theory[END_REF]. The third model assumes the existence of a network of islands and channels which evolves with time and whose tortuosity and connectiveness result in an e ective grain boundary di usivity [START_REF] Raj | Creep in polycrystalline aggregates by matter transport through a liquid phase[END_REF].

Evidence on the nature of the grain boundary structure during the action of pressure solution is di cult to obtain in the laboratory. One of the reasons is, of course, the characteristic time of the process which has to be shortened in the laboratory by changing the chemistry of the uid phase [START_REF] Gratier | Experimental pressure solution-deposition on quartz grains: the crucial e ect of the nature of the uid[END_REF] or by increasing the temperature. For instance, [START_REF] Cox | Experimental dissolution-precipitation creep in quartz aggregates at high temperature[END_REF] have been able to show through an experiment on a ÿne powder of quartz grains at 1200 K that there was indeed a network of islands and channels within the grain contact. Similar results have been obtained at room temperature on rock salt, which should be seen as an analogue material, by [START_REF] Urai | Weakening of rocksalt by water during long-term creep[END_REF] and Spiers et al. (1990). These conclusions have been challenged by [START_REF] Hickman | Kinetics of pressure solution at halite-silica interfaces and intergranular ÿlms[END_REF] who advocated, from their observations of contacts between polished lenses of halite, that no uid phase could be maintained within a contact to allow IPS to occur.

The question of the internal structure of the grain boundary will be solved only by observations collected during the action of IPS. This paper proposes, in parallel to the ongoing experimental e ort [START_REF] Schutjens | Intergranular pressure solution in NaCl: grain-to-grain contact experiments under the optical microscope[END_REF], a numerical modeling of the evolution with time of the grain boundary structure. The solution of initialand boundary-value problems should thus become the mere consequence of accepted assumptions on IPS, which are discussed next, and is used to shed light on the nature of the deformation mechanism.

The ÿrst basic and well-accepted assumption to model IPS is that the solid-uid interaction can be described as a phase transformation between a non-hydrostatically stressed solid and a uid. In addition, the solid is considered here to be isotropic and elastic and the uid in mechanical and chemical equilibrium. The uid phase is composed of a diluent (e.g., water) and a single solute of the same composition as the solid phase. The characteristic length of the problem is less than the grain boundary thickness, so that the di usion time is short compared to the kinetics of the phase transformation. The kinetic law was derived using concepts of irreversible thermodynamics by Lehner and Bataille (1984/85) and then extended by [START_REF] Heidug | Geometrical evolution of stressed and curved solid-uid phase boundaries. Part 1-transformation kinetics[END_REF] to include surface tension e ects. The kinetic driving force is, of course, zero if Gibbs' (1878) equilibrium condition is met.

Little work has been reported so far in the literature on IPS as an initial-and boundary-value problem. The few contributions are concerned with the initial morphology change in time of a solid-uid interface using linear stability arguments. For example, [START_REF] Leroy | Geometrical evolution of stressed and curved solid-uid phase boundaries. Part 2-Stability of cylindrical pores[END_REF] found that the agent driving the instability of a planar surface is the jump in tangential stress component across the solid-uid interface. Surface tension stabilizes the shorter wavelength modes and provides the characteristic length responsible for the selection of a dominant mode of instability. These ÿndings have many similarities with known results on the stability of interfaces driven by surface di usion [START_REF] Srolovitz | On the stability of surfaces of stressed solids[END_REF], which has its roots in the seminal contribution of [START_REF] Mullins | Morphological stability of a particle growing by di usion of heat ow[END_REF]. Stability analyses which di er by the choice of relevant characteristic time should, of course, provide the same neutral stability predictions if they assume the same conditions for thermodynamic equilibrium [START_REF] Grinfeld | Thermodynamics Methods in the Theory of Heterogeneous Systems[END_REF]. Non-linear evolution with time of a solid-uid interface, which is the main objective of this paper, has received less attention apart from the e ort on solidiÿcation and melting processes (e.g., [START_REF] Zabaras | On the calculation of deformations and stresses during axially symmetric solidiÿcation[END_REF]. At the scale of the grain boundary, non-linear evolution due to surface di usion has been considered by [START_REF] Chuang | Non-equilibrium models for di usive cavitation of grain interfaces[END_REF] to capture the relation between time for cavitation and both defect geometry and applied stress. Detailed calculations with Green's functions have revealed the formation of a cusp from planar surfaces or cylindrical cavities [START_REF] Chiu | Stress singularities along a cycloid rough surface[END_REF][START_REF] Yang | Surface morphology evolution in stressed solids: surface di usion controlled crack initiation[END_REF][START_REF] Wang | Shape change of a pore in a stressed solid via surface di usion motivated by surface and elastic energy variation[END_REF].

In the next section, one states the governing equations and their dimensional form and expresses the weak formulation for initial-and boundary-value problems. The driving force for phase transformation is presented once general thermodynamic equilibrium conditions are stated in terms of chemical potential tensors [START_REF] Bowen | Toward a thermodynamics and mechanics of mixtures[END_REF]. The ÿnite-element method is presented next. Its feature is an enhancement of the continuity of the ÿeld variables and the geometrical interpolation on the interface. Further information on the element and its convergence properties are found in Appendix B. The solutions of two initial-and boundary-value problems are presented in Section 4 to shed light on the evolution with time of the grain boundary internal structure. The numerical scheme is validated by comparing the results with the analytical predictions of a linear stability analysis, found in Appendix C.

Governing equations, dimensional analysis and weak formulation

This section pertains ÿrst to the presentation of the physics governing the migration of an interface separating a non-hydrostatically stressed solid and a uid phase in mechanical and chemical equilibrium. The migration results from a phase transformation between a linear elastic, isotropic solid and a uid containing the same material in solution with a diluent. The kinetics of the phase transformation provides the dominant characteristic time of the interface migration which is thus assumed to be large compared to the di usion characteristic time within the uid phase. No di usion occurs within the solid phase. The physics necessary to describe IPS being introduced, a dimensional analysis is proposed prior to the discussion of the weak formulation for initial-and boundary-value problems.

Field equations and kinetic law

The ÿrst objective is to introduce the thermodynamic force responsible for the migration of the solid-uid interface starting from the general thermodynamic equilibrium condition for two phases composed of several constituents. It is shown how this condition includes Gibbs' (1878) equilibrium for the particular solid-uid system of interest. This derivation is found useful to cast our analysis in a general framework and also serves the purpose of deÿning notation.

To deÿne the conditions for thermodynamic equilibrium across an interface, it is convenient to introduce the concept of chemical potential tensor. [START_REF] Bowen | Toward a thermodynamics and mechanics of mixtures[END_REF] was apparently the ÿrst to do so in an attempt to derive a general theory of mixtures without postulating the concept of partial uxes. The deÿnition of [START_REF] Bowen | Di usion in mixtures of elastic materials[END_REF] is amended to express the chemical potential tensor as

K a = a I -A a -1 a ;
(1) in terms of a , a and A a which are the partial speciÿc Helmholtz free energy, the partial mass density and the partial Cauchy stress acting on the ath constituent in a given phase, respectively. Bold letters are reserved for tensorial quantities. The tensor I found in Eq. ( 1) is the second-order identity tensor. The general condition for thermodynamic equilibrium for every constituent present in two phases separated by an interface 9 I , which is discussed by [START_REF] Truskinovsky | The chemical-potential tensor[END_REF], [START_REF] Heidug | Thermodynamics of coherent phase transformation in nonhydrostatically stressed solids[END_REF] and [START_REF] Grinfeld | Thermodynamics Methods in the Theory of Heterogeneous Systems[END_REF], is written here by stating that the jump in normal component of the vector K a • n across that interface must be zero

n • (<K a = • n) = 0: (2)
The bracket stands for the di erence between the function in argument estimated in the two phases across the interface. The vector n is the normal to the stressed interface and its orientation, required in the following, is towards the uid phase.

The deÿnition and the name adopted in Eq. ( 1) are now justiÿed by looking at the particular case of a uid phase. The Helmholtz free energy of a uid phase f is the sum of the partial Helmholtz free energy weighted by the concentration C a of each constituent. The total stress acting locally is the sum of the partial stresses. For a uid in mechanical equilibrium, that partial stress state is hydrostatic and equal to the opposite of the pore uid pressure p f times the ath species concentration C a . In this instance, the chemical potential K a is isotropic and proportional to a + C a p f = a . This scalar is the partial Gibbs free energy of the ath constituent if work can only be performed by volume change. Note that the partial mass density a is nothing but C a f in which f is the mass density of the uid phase. Furthermore, by virtue of Euler's theorem, the partial Gibbs free energy equals the chemical potential of the constituent denoted by a and K a is nothing but a I. These remarks and deÿnitions, classical in mixture theories [START_REF] Bowen | Di usion in mixtures of elastic materials[END_REF], are consistent with Gibbs-Duhem equation in the uid phase:

f + p f f ≡ a C a a + p f f = a C a a ; (3)
and justify the origin of the name attributed to the tensor K a in Eq. ( 1). We now specialize the equilibrium condition (2) to the particular solid-uid system of interest. The uid phase is composed of two constituents or species, the solute and the diluent. Quantities attached to the former and the latter species are marked by a S and a D in subscript, respectively. The solid phase is composed of a single constituent of the same nature as the solute in the uid phase. Field quantities assigned to the solid phase are designated by a s in subscript. The solid-uid interface is envisioned, following the proposition of Gibbs, as a dividing surface of virtual thickness appropriately positioned in the region of rapid variation in densities so that the thermodynamic force associated to the curvature vanishes [START_REF] Gibbs | On the equilibrium of heterogenous substances[END_REF]. This interface could also be seen, however, as a membrane having speciÿc properties and sustaining its own state of stress and deformation [START_REF] Alexander | Thermomechanical equilibrium in solid-uid systems with curved interfaces[END_REF][START_REF] Heidug | A thermodynamic analysis of the conditions of equilibrium at nonhydrostatically stressed and curved phase boundaries[END_REF][START_REF] Heidug | Geometrical evolution of stressed and curved solid-uid phase boundaries. Part 1-transformation kinetics[END_REF] to characterize the energy stored in that region. For the present investigation, the simplifying assumption of an isotropic interface is proposed permitting the mechanical equilibrium conditions across the interface to be written as

t • A s • n = 0; n • A s • n = 2H -p f ;
(4) in terms of the mean curvature 2H , surface tension and two arbitrary vectors t tangent to the interface. Eqs. (4) are the classical Laplace conditions which rely on the assumptions that the uid is in equilibrium (A f = -p f I) and thus disregard any shear e ect along the interface due, for example, to the uid viscosity. Combining Eqs. ( 1) and ( 4) and the deÿnition provided above for the solute chemical potential with the jump condition (2) provides s +

p f -2H s -S = 0;
(5) which is Gibbs' famous thermodynamic equilibrium condition [START_REF] Gibbs | On the equilibrium of heterogenous substances[END_REF] for the solute in contact with the solid phase.

Having deÿned the concept of thermodynamic equilibrium, we now turn our attention to the class of kinetic laws which govern the migration of the solid-uid interface under conditions which remain close to the equilibrium condition (5) to render applicable the results of irreversible thermodynamics. This discussion is restricted to our particular solid-uid problem. The solid phase occupies at time t the domain (t) called the current conÿguration. If the solid phase is released from any external loading, its position then deÿnes the reference conÿguration denoted 0 (t). A point X in the domain 0 (t) is mapped to x in (t) by a ÿnite transformation. Note that dissolution and precipitation change the total mass contained in 0 (t) which is thus also time dependent. This change in mass results from the mass ux J n through the solid-uid interface 9 I (t) and is deÿned per unit area of the interface in its current conÿguration by

J n = s (ẋ • n -S n ); (6)
in terms of S n and ẋ which are, respectively, the speed of displacement of the soliduid interface and the material velocity of the point x currently on the interface. It is more convenient in our analysis to employ the mass ux J N per unit area of the interface in the reference conÿguration. To relate the two uxes, let one consider a patch on this interface of areas dS and dS 0 and oriented by the normals n and N over the current and reference conÿgurations, respectively, and take note that J n dS = J N dS 0 ≡ -0s S N dS 0 :

(7)

The deÿnition of J N relies on the Lagrangian speed of propagation S N and the reference value of the solid density 0s (Truesdell, 1977, Chapter II.6). According to Eq. ( 7), the ux within the uid phase is positive if the non-material velocity S N is negative signalling a dissolution process. A precipitation requires a positive velocity S N and corresponds to an interface migrating towards the uid phase with a negative ux J N .

The thermodynamic driving force conjugate to the mass ux J n denoted by is shown by Lehner and Bataille (1984/85) and [START_REF] Heidug | Geometrical evolution of stressed and curved solid-uid phase boundaries. Part 1-transformation kinetics[END_REF], from an estimate of the dissipation during migration, to be the left-hand side of Eq. ( 5), which is the expression taken by the jump condition (2) for the solute and the solid phase in our problem. Any function L( ) having the following property:

J n = L( ) with L( ) ≥ 0 and = s + p f -2H s -S ; (8)
is a candidate for a kinetic law.

To determine the non-material velocity, or speed of propagation, S N of the interface with Eqs. ( 8) and ( 7) one needs to know the solute speciÿc chemical potential. It is a function of the pore uid pressure and the solute concentration. In [START_REF] Ghoussoub | Solid-uid phase transformation within grain boundaries during compaction by pressure solution[END_REF], this function is derived for a mixture of ideal solutions and the ÿnal expression reads

S (p f ; C S ) = * S + p f -p 0 f S + RT M S log C S C 0 S ; (9)
in which * S should be seen as a reference chemical potential taken at a pore uid pressure p 0 f and a solute concentration C 0 S . The constants R, T and M S are the universal gas constant, the temperature and the solute molar mass. The problems treated in what follows are partitioned into two classes. The ÿrst class includes open systems: the uid phase exchanges mass with a remote reservoir of inÿnite dimension so that the uid pressure and the solute concentration are una ected by local changes in morphology of the solid-uid interface. These two quantities are constant during the whole analysis as well as, consequently, the solute chemical potential according to Eq. ( 9). In the second class of problems, the systems are closed: the uid phase is trapped within the grain boundary. In such closed systems, the diluent mass is preserved and the pressure and the solute concentration vary with the local evolution of the solid-uid interface. In the course of this evolution the current volume of the uid phase V f and the mass of the solute m S vary. It involves solute di usion within an isolated cavity at a rate su ciently fast for the uid phase to be still considered homogeneous at all times. Assuming an ideal solution and a constant isothermal compressibility Ä D for the diluent, it is shown in Appendix A that the volume of the uid phase is the following function of the solute mass and uid pressure:

V f (p f ; m S ) = m S p f p 0 f V 0 S m 0 S + V 0 D exp(-(p f -p 0 f )Ä D ); (10)
in which the superscript 0 identiÿes the reference state and V D denotes the volume of the diluent. In our analysis, the volume of the uid phase as well as the solute mass are known and Eq. ( 10) provides the means to compute the pore uid pressure. The pore pressure and the solute concentration (C S = m S =(m S + m 0 D )) determine the solute chemical potential (9) which is required to compute the speed of propagation pointwise along the solid-uid interface.

Further assumptions and dimensional analysis

The objective of this second subsection is to provide the governing equations in a dimensionless form with the assumption of inÿnitesimal transformation (inÿnitesimal strain and displacement). While doing so, further assumptions are provided to complete the deÿnition of our problem.

The solid phase is assumed to be free of stress in its reference conÿguration and to be homogeneous, isotropic and linear elastic. Any new material deposited on the solid phase from the uid phase is assumed to be added coherently (no build up of residual stresses) and to have the same elastic properties as the material initially in a solid state. The elastic properties of the solid phase are thus described by G and K, the modulus of elasticity in shear and bulk deformation, respectively. A superposed tilde over any constant or ÿeld variable signals that this physical quantity has dimension. This convention was, of course, not adopted in the previous part of the section. Every dimensionless variable of generic name A(x; t) is obtained by dividing the physical quantity Ã(x; t) by the reference value A R . Reference values will henceforth be denoted by a subscript letter R. These deÿnitions, which are summarized in Table 1, are now discussed.

Table 1 Characteristic or reference values for physical quantities required for the dimensional analysis, assuming a linear kinetic law as in Lehner and Bataille (1984/85) for a silica-water reaction at a temperature of 217

s =(L R R ) over L R under a driving force X R Mass ux J R Deÿned as R L R =t R 414 kg=(m 2 s) Kinetics constant L R Deÿned as 2 R L R =(t R R ) 1 :1 × 10 -14 kg s=m 4
The mass density reference value R is set to the solid mass density in its reference conÿguration 0s for the particular case of quartz. Any stress à is normalized by the characteristic stress R which is chosen as the magnitude of the lithostatic pressure in the Earth's crust at 1 km depth. Surface tension ˜ has dimension of stress times length and the reference R , used by [START_REF] Leroy | Geometrical evolution of stressed and curved solid-uid phase boundaries. Part 2-Stability of cylindrical pores[END_REF], is 0:1 Pa m. Dividing this quantity by the reference stress provides the characteristic length L R of 10 -8 m.

The only physical time scale discussed so far enters the kinetic law. Here, we select the simplest kinetic law of the class (8) corresponding to a linear relation between the mass ux and the thermodynamic force with a phenomenological parameter L. From the data of [START_REF] Rimstidt | The kinetics of silica-water reactions[END_REF], the value 1:1×10 -14 kg s=m 4 is assigned to this parameter, at a temperature of 217

• C. This value is the characteristic L R which is used to deÿne the characteristic time in Table 1 (equal approximately to a week). The dimensional analysis results in the following relation between dimensionless non-material velocity and thermodynamic force

S N = -: (11)
This result is obtained with the assumption of small transformation which stipulates that the reference and current conÿgurations are su ciently close to disregard to ÿrst order the di erence between dS and dS 0 in Eq. ( 7). The dimensionless driving force introduced in Eq. ( 8) with the small transformation assumption reads

= 1 2 K -2 3 G (tr U) 2 + G tr(U 2 ) + (1 + tr U)(p f -2H ) -S ; (12)
in which U is the inÿnitesimal strain tensor. Note that Eq. ( 12) was obtained by replacing the ratio ( ˜ = R ) by its value of one. Note also that the Helmholtz free energy of the solid phase under isothermal conditions is simply the elastic stored energy which is of same order as the variation in the mass density found in the numerator of the second term on the right-hand side of Eq. ( 8). This comparison justiÿes our choice to keep the inverse of the change in volume in Eq. ( 12) with the small transformation approximation. The importance of retaining this term in a linear stability analysis was discussed by [START_REF] Leroy | Geometrical evolution of stressed and curved solid-uid phase boundaries. Part 2-Stability of cylindrical pores[END_REF] and was already clear from the work of Lehner and Bataille (1984/85).

Weak formulation

The last topic of this section concerns the presentation of the weak formulation of the initial-and boundary-value problem. The elastic material occupies the volume (t) with boundary 9 (t) and it is in contact with a uid phase over the interface 9 I (t). The phase transformation acts on this interface 9 I (t) which evolves with time and its position constitutes one of the unknowns of the problem. Note that the uid phase is in chemical and mechanical equilibrium and thus always homogeneous for both open and closed systems. For this reason, a numerical discretization is only required for the solid phase.

At any time t at which the shape of the solid-phase reference conÿguration is known, the principle of virtual work provides the mechanical equilibrium condition

(t) (K -2 3 G)tr(U)∇ • u + 2GU : ∇ u dV + u • N(2H -p f ) dS = 0; (13)
in which u and ∇ are a kinematically admissible virtual displacement and the gradient operator, respectively. Note that the applied traction vector over the solid-uid interface accounts for the e ects of surface tension according to Laplace's conditions (4). The solution of Eq. ( 13) provides, for a known position of the solid-uid interface, the material displacement u over the physical domain. This information is then used to determine the Helmholtz free energy of the solid at the interface and thus to estimate the thermodynamic driving force (12) responsible for its migration. The non-material velocity of the interface is then provided by the kinetic law (11), which cannot be enforced pointwise if a numerical approximation is introduced. The error or residual is (S h N + ) in which S h N is the numerical approximation of the speed of propagation. We impose here that the weighted average of this residual integrated over the interface vanishes

9 I (t) S h N + 1 2 (K -2 3 G)(tr U) 2 + G tr(U 2 ) + (1 + tr U)(p f -2H ) -S ' dS = 0 (14)
for some appropriate choice of the weight function '. The solution of Eqs. ( 13) and ( 14) provides the two unknowns of the problem which are the displacement over the whole domain and the non-material velocity S h N of the solid-uid interface. The numerical scheme designed for solving the evolution problem and based on this weak formulation is found in the next section.

A class of ÿnite elements for IPS

The need for a special ÿnite-element scheme to capture the evolution with time of the solid-uid interface can be motivated from the interpretation of numerous linear stability analyses found in the literature [START_REF] Srolovitz | On the stability of surfaces of stressed solids[END_REF][START_REF] Leroy | Geometrical evolution of stressed and curved solid-uid phase boundaries. Part 2-Stability of cylindrical pores[END_REF]. The results presented in Appendix C on the stability of a rectangular block point to the same conclusion. An interface between a stressed solid and a uid in mechanical equilibrium is destabilized if the stress component in the direction tangent to the interface di ers from the uid pressure. The stabilizing factor, due to surface tension, is proportional to the mean curvature. The way this curvature is estimated over the discretized interface is now discussed. A classical ÿnite-element scheme relies on a C 0 interpolation rendering the tangent vector along the solid-uid interface discontinuous across elements. Any stress concentration close to the element boundaries is then destabilizing and cannot be balanced by surface tension if the measure of curvature does not account for the discontinuity at element boundaries. This problem was found to be alleviated by a C 1 interpolation of the solid-uid interface. A rapid local variation of the tangent a ects the curvature of the elements in the neighborhood and thus triggers the stabilizing in uence of surface tension. A similar idea was used by [START_REF] Yang | Surface morphology evolution in stressed solids: surface di usion controlled crack initiation[END_REF] in a boundary integral scheme.

This section is made of three parts. The ÿrst describes the class of elements considered with a C 1 interpolation along the solid-uid interface. The second deals with the time discretization and the solution algorithm. The basic unknowns are the material displacements (and its surface gradient) and the position (including the tangent) of the interface. The third part concerns the mesh adaptation which is used to parameterize the interface as well as to modify the mesh over the reference conÿguration whose boundary evolves with time. Proof of convergence of the numerical scheme is discussed in Appendix B. The reader most interested by the physics of the problem can skip this section altogether. Note that although the governing equations were deÿned in the previous section in 3 dimensions, the ÿnite-element scheme to be presented is restricted to the plane. Extension to 3D while straightforward from the theory remains a major technical challenge.

Elements with C 1 interpolation on part of their boundary

Our objective is the construction of a ÿnite element having C 1 continuity on one of its sides, if located at the solid-uid interface, and having the usual C 0 properties in the bulk of the solid phase. The starting point is the interpolation along the interface which is divided in a number of elements. Any point X on that interface is located on the boundary of an element (Fig. 1) and has its position interpolated by

X( 1 ) = 4 a=1 H a ( 1 ) Ỹa ; (15)
in which Ỹa are generalized nodal position vectors comprising 2 nodal positions (Y 1 ; Y 2 ) and two tangent vectors ( 1 T 1 ; 2 T 2 ), having for norm 1 and 2 . These functions H a ( 1 ) are ÿrst-order Hermite polynomials of the natural coordinate 1 . The scalars 1 and 2 are the Jacobian at 1 equal to -1 and +1 of the transformation between the curve (s) of arc length s and the ( 1 ) axis (Fig. 1a). These scalars are part of the interface parameterization and it is proposed to determine them from constraints of the type = G ( Ỹa ; Ỹ3 ):

The two unknowns in Eq. ( 16) are deÿned in terms of an extra node Y 3 , shown as an open circle in Fig. 1a. For example, the function G could be constructed such that Y 3 be positioned on the interface at the coordinate 1 equal zero. A discussion of this particular choice is found in [START_REF] Ghoussoub | Solid-uid phase transformation within grain boundaries during compaction by pressure solution[END_REF].

Interpolation ( 15) should be compatible with the ÿnite-element discretization over the solid phase having part of its boundary on that interface (Fig. 1b). For that purpose, a special class of elements is constructed having this speciÿc property. The elements have their shape functions constructed by combining Lagrange's and Hermite's polynomials. Their exact expressions are found in Appendix B and it su ces here to know that the interpolation over the element is given by

X( 1 ; 2 ) = N a=1 N a ( 1 ; 2 ) Xa ; (17)
in which N a and Xa are the shape functions and the generalized nodal coordinates (Fig. 1b). The list of generalized coordinates is (X 1 ; X 2 ;Y 2 ; Y 1 ; 2 T 2 ; 1 T 1 ; X 7 ; X 8 ; X 9 ; X 10 ) if the number of nodes N is set to 10 and according to the connectivity deÿned in Fig. 1c. The next step in the construction of the ÿnite-element scheme concerns the interpolation of the displacement ÿeld. Here, we favor the idea of an isoparametric element in which displacement and position interpolation are thus identical. The collection of generalized nodal displacement ũa includes the material displacement as well as its surface gradient: (u 1 ; u 2 ; u 3 ; u 4 ; 2 (∇ u) 2 ; 1 (∇ u) 1 ) for the basic 4 noded element and with ∇ deÿning the surface gradient over the interface . In view of the deÿnition of the scalars ( = 1; 2) as Jacobians of the transformation between (s) and ( 1) at 1 equal to -1 and +1, the surface gradient (∇ u) times could be replaced by (∇ u) , the gradient of the displacement over the -axis.

The mesh adaptation scheme

The evolution in time of the solid-uid interface requires the modiÿcation of the ÿnite-element discretization of the domain occupied by the solid phase. The solution adopted to deÿne the mesh over the reference conÿguration 0 (t) relies on a boundary conforming mapping. It is a variational method in which the boundaries of the reference conÿguration, the physical domain, coincide with the transformed coordinate lines of the computational domain which are equi-potential. This method, based on the work of [START_REF] Brackbill | Adaptive zoning for singular problems in two dimensions[END_REF], aims to optimize three properties of the discretization which are its smoothness, the orthogonality of the coordinate lines and the variation or concentration of the elements in regions of interest. The -coordinate line mapped on the computational domain has a gradient deÿned by the vector 9 =9X i . Smoothness of the coordinate lines spacing is then measured by the sum of the squares of the norm of the gradient vector for the two lines

I S = 0 9 9X j 9 9X j dV: (18)
Orthogonality is estimated by the square of the scalar product of the gradient between the 1 -and 2 -lines

I O = 0 9 1 9X j 9 2 9X j 2 J 3 dV; ( 19
)
in which J is the Jacobian of the transformation from the computational to the physical domain. The concentration of the discretization is measured by

I C = 0 !J dV; (20)
with the introduction of !, a function of the Helmholtz free energy [START_REF] Ghoussoub | Solid-uid phase transformation within grain boundaries during compaction by pressure solution[END_REF]. The Euler-Lagrange equations of the adaptation scheme are then obtained by minimizing the linear combination of the 3 potentials in Eqs. ( 18) -( 20)

I S + O I O + C I C ; (21)
where O and C are two scalars. These two scalars could be seen as Lagrange multipliers since I O and I C do not lead to elliptic operators at a di erence with I S [START_REF] Brackbill | Adaptive zoning for singular problems in two dimensions[END_REF]. Furthermore, the selection of these scalars is dictated by the fact that the 3 potentials in Eqs. ( 18) -( 20) do not share the same dimension, as pointed out by [START_REF] Kreis | Application of variational method for generating adaptive grids[END_REF]. If I S is of order 1, then I O and I C are of order h and !h, respectively, in terms of the parameter h which is the characteristic element size. Appropriate scaling in Eq. ( 21) requires O and C to be of order h -4 and ! -1 h -4 , respectively.

Various di culties arise with the use of this class of variational methods for adapting numerical discretization which are discussed in [START_REF] Christodoulou | Discretization of free surface ows and other moving boundary problems[END_REF] and [START_REF] Tsiveriotis | Boundary-conforming mapping applied to computations of highly deformed solidiÿcation interfaces[END_REF]. The latter contribution is closer to our problem since it was concerned with the migration of solidiÿcation fronts. Among those numerical problems, we should mention the di culty to concentrate -coordinate lines parallel to convex boundaries compared to concave boundaries and the absence of parameterization of the curved solid-uid interface. Solutions to these problems and further information on the adaptation scheme, including the linearized version of the Euler-Lagrange equations for the class of ÿnite elements proposed above, are found in [START_REF] Ghoussoub | Solid-uid phase transformation within grain boundaries during compaction by pressure solution[END_REF]. The values chosen for every parameter of the mesh adaptation in the calculations presented in this paper are also found in that reference.

Time discretization and solution strategy

Time is discretized and it is assumed that at time t n the reference conÿguration O (t n ) is known. The four steps now presented summarize the procedure proposed to compute the next position of the solid-uid interface at time t n+1 .

First, knowing the reference conÿguration at time t n one determines the material displacement through the application of the principle of virtual work (13). The strain and stress are then known over the domain O (t n ) and, in particular, over the solid interface 9 I (t n ) where the thermodynamic force has to be computed. In a second step, the weak form of the kinetic law ( 14) is solved to provide the non-material velocity S h N (t n ). The position of the solid-uid interface is updated in a third step by an explicit scheme based on the non-material velocities just computed. It remains in the fourth step to modify the parameterization of the interface and, subsequently, to adapt the ÿnite-element discretization over the solid phase. These two mesh adaptations are not conducted at every time step and a simple criterion is proposed to activate this last part of the procedure: the mesh is adapted over the volume or the interface only if there is at least one element whose area increases by more than a certain percentage, set to 3, for most of the results to be presented. The accuracy of the explicit scheme has been tested by comparing several choices of time step for the same problem. It has been found that the solution is accurate if during a time step any element does not have a relative change in its surface exceeding 10 -2 %. It is clear that an implicit scheme could be more interesting since it would allow larger time step. It would require the simultaneous resolution of the mesh adaptation scheme, the weak form of the kinetic law and the principle of virtual work. However, the implicit scheme does not have a symmetric tangent operator in the absence of a variational formulation for the coupled problem. Furthermore, the increase in number of nodal unknowns compared to our staggered algorithm may cancel the beneÿt of an implicit scheme.

The ÿnal word on the numerical scheme concerns the case of a constrained uid phase for a closed system. In that instance, the solution procedure has to be amended to permit as a ÿfth step the computation of the chemical potential at time t n+1 . For that purpose, the new volume of the inclusion is ÿrst computed as well as the total mass present in the solid phase. The di erence in mass with respect to the previous time step determines the variation in solute mass in the uid phase and thus the current concentration C S . The pore uid pressure is then determined by solving the non-linear equation (10). Note that the time step is always smaller for the study of closed systems in order to avoid spurious uctuations in pore pressure.

Results

The main objective of this paper is to gain a better understanding of the evolution of the grain boundary structure during a loading cycle. In Fig. 2, a transmission optical micrograph of a grain boundary in rock salt after deformation by IPS in the laboratory is presented (result obtained in the HPT Laboratory, Faculty of Earth Sciences, Utrecht University). The grain boundary shows the typical island and channel structure discussed before. The channels are tubular cavities, about 3 m in diameter, which were uid ÿlled during the experiments. Note also the existence of isolated inclusions within the grain boundary. The contact made by the islands are likely to be dry regions even though the existence of thin ÿlms, especially during the experiments, cannot be ruled out. This experimental observation motivates the proposition of the idealized grain boundary depicted in Fig. 3. Two features in this idealized grain boundary are sketched. They will be at the basis of two initial-and boundary-value problems proposed to model the time evolution of an idealized grain boundary.

The ÿrst structure considered presupposes that a thin uid ÿlm penetrates at least partly the grain boundary and provides an e cient channel for the transfer of any dissolved material to the pore uid. The corresponding model problem is referred to as problem 1 in Fig. 3 and consists of a rectangular block in contact on its top side with the uid phase. Only the case of a constant solute chemical potential is explored, corresponding to the deÿnition of an open system in Section 2.

The second microstructural feature found in Fig. 3 is a cylindrical uid-ÿlled tube trapped within the grain boundary. This geometry is of course reminiscent of the channel structure observed in Fig. 2. The model problem extracted from that microstructure and labeled problem 2 in Fig. 3 consists of a cylindrical cavity in an elastic solid sustaining a remote compression normal to the grain boundary. The energy associated with the solid-solid contact is disregarded and a circular cross section is the equilibrium shape. For reasons of symmetry, only a quarter of the structure is studied. To discuss the in uence of the uid phase, both the open and closed systems are considered. For open systems, the solute chemical potential is kept constant and takes the value necessary to have thermodynamic equilibrium under the initially hydrostatic stress conditions. For closed systems, the uid phase is assumed isolated within the grain boundary and the chemical potential varies according to the uid pressure and the solute concentration calculated in the plane of study.

The rectangular block

The ÿrst problem considered consists of the rectangular block, thickness H and length L, shown in Fig. 4. The solid-uid interface BB is initially plane and the compression of the block is due to a force F or a displacement U prescribed on the right lateral boundary A B . The force is imposed by a rigid plate of constant length denoted H G . The thickness H is always less than H G and the uid pressure acts on the top of the plate over the length H G -H . This particular choice was found convenient for the linear stability analysis. Conditions of plane strain deformation prevail in the elastic block. The numerical data are summarized in Table 2. The fundamental solution is ÿrst discussed prior to the analysis of the initial morphology evolution of the solid-uid interface when its position is perturbed by a cosine wave of small amplitude. The long term evolution of this type of perturbation is then considered.

The fundamental solution of this problem corresponds to an initially rectangular block whose shape remains homothetic at all times. The rate of change of the thickness H (t) is equal to the non-material velocity S N of the solid-uid interface, which, according to the linear kinetic law (11), reads

Ḣ (t) = - 1 - 4G ( 0 + p f ) 2 -p f + S ; (22)
in which and 0 are Poisson's ratio and the uniform stress acting in the x-direction, respectively. This stress 0 depends on the type of boundary conditions

according to 0 = 2GU 1 - -p f 1 - for displ: control; 0 = -p f + F + p f H G H for force control: (23)
The stress 0 is thus independent of H for displacement control: the block thickness either increases or decreases at a constant rate, as illustrated by the two dotted curves in the phase diagram of Fig. 5, for the particular choice of 0 equal to -p f . A single value of 0 permitting thermodynamic equilibrium is found from Eqs. ( 22) and ( 23), if the pore uid pressure is smaller than the dimensionless solute chemical potential.

If that last condition is not met, Ḣ is always negative and dissolution occurs at a constant rate. If p f is smaller than S , dissolution or precipitation occurs depending on the e ective stress ( 0 + p f ) found in Eq. ( 22). The importance of this e ective stress ( 0 + p f ) will be conÿrmed by the stability analysis. For force control, the existence of an equilibrium state, marked by a solid dot in Fig. 5, also requires the pore pressure to be less than the dimensionless chemical potential of the solute S . The importance of the order relation between p f and S is illustrated by two solid curves for force control in Fig. 5. If the pore pressure is smaller than S , the thermodynamic driving force does change sign and an equilibrium is reached for a thickness which is used as a reference value H Ref . Fig. 5 is also instrumental in explaining the linear stability of an equilibrium following a homothetic perturbation of the block which remains thus rectangular in shape. For displacement control, this equilibrium is neutrally stable since the fundamental solution is independent of H . For force control, the equilibrium is unstable since any inÿnitesimal increase in the block thickness leads to a positive value of Ḣ which results in further increase in the block thickness. On the contrary, a perturbation which would reduce the block thickness renders the rate Ḣ negative and results in further reduction of the grain thickness, away from the equilibrium value. The rate of dissolution or precipitation is independent of H in a displacement control test (DC), as shown by the dotted lines for the special case of 0 = -p f . For the force control case (solid curves), an equilibrium (solid dot) to exist requires the dimensionless chemical potential to be larger than the pore pressure. Note that this equilibrium is unstable for a homothetic perturbation.

4.1.1. Initial evolution and validation of the numerical scheme One of the reasons for studying problem 1 is that the initial evolution of wavy perturbations, of small amplitude compared to the block size, can be compared to the prediction of a linear stability analysis. This linear stability analysis, summarized in Appendix C, provides the following expression for the stability exponent:

= -2 + ( 0 + p f ) 2 (1 -) G cosh 2 ( H ) cosh( H ) sinh( H ) + H ; ( 24
)
in which is the wavenumber of the perturbation, always a multiple of =L. Note that the same result (24) applies for force and displacement control. If this stability exponent is positive for a given the solid-uid interface is unstable. Conversely, a negative for all deÿnes a stable interface. The analysis of this stability exponent is started with the ÿrst terms -2 on the right-hand side of Eq. ( 24) which is due to surface tension. It is negative and thus has a stabilizing e ect for su ciently large as can be seen from the four curves drawn in Fig. 6, for two choices of H and of the compressive stress 0 . The second term on the right-hand side of Eq. ( 24) is always positive or equal to zero and is thus destabilizing. It involves the product of and three factors. The ÿrst factor is Fig. 6. Stability exponent as a function of the normalized wavenumber for two values of the thickness H and of the compressive stress 0 . The stability exponent is always negative, and the interface stable, if 0 is equal to minus the uid pressure. If this condition is not met, there is a range of wavenumbers for which the stability exponent is positive and the interface is thus unstable. In the inset, the stability exponent is drawn in the long-wavelength limit for comparison with the analytical result obtained for = 0 (homothetic mode).

proportional to the square of the e ective stress ( 0 + p f ) which needs to di er from zero for the interface to be unstable. As in [START_REF] Leroy | Geometrical evolution of stressed and curved solid-uid phase boundaries. Part 2-Stability of cylindrical pores[END_REF], it is the di erence between the interface-parallel compressive stress and the pore pressure which controls instability. If 0 is set to -p f , the stability exponent is always negative, as seen in Fig. 6 (curve 4). There is thus for any 0 di erent from -p f a range of for which the stability exponent is positive because surface tension e ects are then of second order in and not su cient to stabilize the interface. The wavenumber for which the stability exponent is maximum is called the dominant wavenumber. The present discussion has been conducted so far disregarding the in uence of the third factor in the second term on the right-hand side of Eq. ( 24), which is due to the ÿnite size of the studied domain and is a function of the normalized wavenumber H . This geometrical function is close to one for values of H su ciently large and its in uence can then be disregarded. This can be seen from Fig. 6 by comparing the dominant normalized wavenumber for the same value of 0 = -5p f but di erent block sizes (curves 1 and 3). According to Eq. ( 24), if the geometrical function of H is set to 1, the dominant is then proportional to the square of the e ective stress. The dominant dimensionless wavenumber in Fig. 6 should then scale by H . This is the case here with an error margin of 5%, justifying the statement above on the in uence of the ÿnite geometry. This geometrical function of H has, however, a second salient feature: it is singular for tending to zero leading to a ÿnite value of in that limit, as illustrated in the inset of Fig. 6. This limit is always positive if 0 di ers from -p f and inversely proportional to H . It is smaller than the maximum value of conÿrming that the wavy perturbations are dominant, at least for the choice of 0 made in Fig. 6. A word of caution concerning this limit is now in order since the long-wavelength limit for is found to be independent of the type of boundary conditions speciÿed. This is at variance with the ÿndings on the stability of the fundamental solution for a homothetic mode ( = 0) discussed earlier. The reason for that discrepancy, discussed in Appendix C, is a singular limit in the boundary conditions of the stability analysis.

We now turn our attention to the comparison between the numerical solution and the linear stability prediction of the evolution of the solid-uid interface perturbed by a cosine with wavenumber =L and amplitude 5% of H . The block has a square shape with H set to either 500 or 1000. The displacement is prescribed on the side resulting in a compressive stress 0 , deÿned in Eq. ( 23), of either -1 or -5 times p f for a perfect block. Results are presented in Fig. 7 in which the normalized distance function d 1 (surface between the actual position of the interface and the perfect interface at equilibrium) is plotted as a function of time (solid curves). For 0 = -p f , the distance function is decreasing with time due to an evolution of the interface towards the at equilibrium shape. To the contrary, for 0 = -5p f , the distance increases signaling the development of the instability. Note that in all cases, the rate of growth or decay of the perturbation coincides with the prediction of the linear stability analysis (dotted Fig. 8. The non-material velocity at the point along the interface having largest curvature as a function of time for the square block under compression presented in Fig. 4 and having a length of 100. In (a), the horizontal stress 0 is -18 and the long-wavelength imperfection has an amplitude of 1%, 2% and 5% times the block length. In (b), this amplitude is set to 1% and the stress is set between -18 and -22. The reference stress is 10 MPa. A time of 1000 and a velocity of 1 represent 20.3 years and 0.5 m=yr, respectively. lines) which are also found in Fig. 7. This agreement between numerical and analytical results validates the ÿnite-element scheme proposed in the previous section.

The long term evolution: formation of a dissolution pit

The non-linear evolution of the solid-uid interface of problem 1 is conducted ÿrst for a square block of length 100 and a displacement control resulting in a compressive stress 0 deÿned in Eq. ( 23) of -18. The results of the computation are presented in Fig. 8a where the non-material velocity at the point of maximum curvature (right-hand side of the interface in Fig. 4) is a function of time for an amplitude of the imperfection set to 1%, 2% and 5% the block length. The main observation from this plot is the existence of two stages during the history of the velocity. The evolution is ÿrst rather slow during a period which lasts up to 4000 times the characteristic time. Then, the dissolution (negative velocity) becomes important and the non-material velocity reaches a value which is two orders of magnitude larger than in the ÿrst stage. The duration of stage 1 is a ected by the imperfection size (Fig. 8a), as well as the value of the stress 0 (Fig. 8b). It is seen that an increase in the magnitude of 0 by 20% results in a reduction of the ÿrst stage duration by nearly an order of magnitude.

Similar results are presented as solid curves in Fig. 9 for a square block of length 500, a stress 0 of -18 and an imperfection amplitude of 1%, 2% and 5% the block length. The transition between the two stages is sharper and, as before, the increase in the non-material velocity magnitude does not appear to become bounded. To understand this drastic evolution of the dissolution process, the position of the interface has been Fig. 9. The evolution with time of the non-material velocity at the point along the interface having largest curvature, for the rectangular block of Fig. 4 having a length of 500, is presented for various imperfection sizes, a stress 0 of -15 (solid curves) and for a stress of -18 and an imperfection size set to 5% times the block length (dashed curve). The dotted curve is obtained for conditions similar to the ones selected for the dashed curve except for the chemical potential which is larger, resulting ÿrst in a deposition over the whole interface. A time of 1000 and a velocity of 1 represent 20.3 years and 0:5 m=yr, respectively. plotted in Fig. 10 at 17 times identiÿed by the letter T in Fig. 9 for the particular case of an imperfection of 1%. At time T 0 , the interface has the cosine shape of the imperfection. Between time T 0 and time T 3 , the morphology evolution of the interface is rather di use. The amplitude of the imperfection increases with deposition on the left and dissolution on the right. This is due to the variation in thermodynamic force along the interface resulting from a larger Helmholtz energy in the thinner section of the block. Note that between T 0 and T 3 there is a shift to the right of the in ection point along the interface which is where, initially, the thermodynamic equilibrium condition is met. This shift is due to an acceleration of the dissolution process compared to the deposition process. This acceleration is more obvious between time T 3 and T 6 , just before the transition between stages 1 and 2 as deÿned above (Fig. 10). At time T 6 , the morphology evolution is marked by a rapid dissolution over 10% of the block length. The rate of dissolution remains, however, of the same order as the rate of deposition. The important increase in dissolution occurs only after that time and takes place over a section which decreases in size, as seen in the inset of Fig. 10. Note that at time 16, the dissolution is now taking place over a length twice the characteristic length which is due to surface tension. However, the rate of dissolution keeps on increasing Fig. 10. The position of the interface at 17 times deÿned in Fig. 9 with the letter T. The zoom on the dissolution pit in the inset reveals that, at time 16, the dissolution is localized on a length scale which is of the same order as the characteristic length due to surface tension (10 -8 m). and surface tension does not have the stabilizing in uence found in the linear analysis. No physical mechanisms are accounted for here which could stop the formation of this sharp dissolution pit or cusp.

The reason for this unbounded growth can be explained by analyzing Fig. 11 which provides two types of information in the region around the tip of the dissolution pit. On the right-hand side of this ÿgure, the ÿnite-element mesh over the reference conÿguration at the end of the simulation is presented for the case of a block with length 100, a stress 0 of -18 and a 2% imperfection. The L at the center provides the scale and the horizontal and vertical bars have for length 1:8 and 0:9 (×10 -8 m), respectively. These numbers allow to estimate the curvature at the tip of the cusp as well as the size of the smallest element. Note that this element has changed its size by a ratio of approximately 40, thanks to the mesh adaptation scheme. On the left-hand size of Fig. 11, 26 isocontours of equivalent Von Mises stress are presented in the range of 15-50 (×10 MPa). In the zone close to the tip covered by nine elements, where isocontours are not plotted, this dimensionless equivalent stress reaches the value of 150 (×10 MPa). The shape of those isocontours is reminiscent of the stress distribution at the tip of a sharp crack. They illustrate the importance of the stress concentration at the tip of the dissolution pit and the role played by the Helmholtz free energy. It is that part of the driving force which now controls the rate of dissolution. The simulations have been stopped when the formation of the dissolution pit was clear and no physical argument was proposed to assess the local failure of the material. One could think of a fracture mechanics argument based on an estimate of the stress intensity factor or invoke a plasticity yield criterion, in view of the large stresses mentioned above. For quartz tested in the lab, fracture seems to be the observed failure mechanism. However, creep plasticity is often invoked to be active in nature due to the presence of minute amount of uid within the solid phase. Here, it su ces to say that IPS has activated a second mechanism which is contributing to irreversible modiÿcations of the local morphology which are not captured by the present scheme.

The physical interpretation of the sudden acceleration in the dissolution rate seen in Figs. 8 and 9 is now clear. It is a localization of the phase transformation in a small portion of the block resulting in a dissolution pit or cusp. The latency period, deÿned by the duration of stage 1 deÿned above, is function of the imperfection size and the stress 0 .

The role of the chemical potential on the local morphology change is now discussed. It is a constant scalar selected to have thermodynamic equilibrium pointwise along the interface for the perfect block. This value has been changed to 3.032, corresponding to an increase of 2=3 of a percent with respect to the equilibrium value, to trigger a general deposition over the whole interface. The block has a length of 500, an Fig. 12. Evolution of the solid-uid interface for the square block presented in Fig. 4 and having a length of 500 and an imperfection amplitude 5% this length. The compressive stress 0 is -18 (×10 MPa). The chemical potential is selected to induce initially an overall deposition. However, after time t 7 , the velocity at the point of largest curvature is reversed signaling dissolution through the newly added material. The various times are deÿned in Fig. 9 with the letter t.

imperfection of 5% and the stress 0 is -18. The position of the interface at 15 times are presented in Fig. 12. These times are identiÿed by the letter t on the dotted curve of Fig. 9. The general evolution of the interface is rather similar to the one described above except for a single feature: the amplitude of the imperfection is increasing with the addition of a translation due to the general deposition. The increase in amplitude is due to the build up of a stress concentration in the thinner section of the block. This concentration is su cient at time t 7 to balance the extra contribution of the chemical potential to the driving force and to stop the deposition process on the right-hand side of the block. From that time on, the dissolution accelerates while deposition continues on the left-hand side. This dissolution occurs through the newly added material until time t 11 . After that time, the material initially present is also dissolved. Note from Fig. 9 that t 11 is already in stage 2 corresponding to the localization process. Note also from that ÿgure that the latency period before localization is little a ected by the change in chemical potential by comparing the dotted and dashed curves, the latter being obtained with the chemical potential necessary for equilibrium of a perfect block. The conclusion of this simulation is that a supersaturated pore uid does not rule out local dissolution processes driven by stress concentrations.

The cylindrical tube as an open system

The next two subsections deal with model problem 2 for a cylindrical uid-ÿlled cavity subjected to the loading shown in Fig. 13. The energy associated with the solid-solid interface (along the horizontal axis of symmetry) is disregarded and the thermodynamic equilibrium shape of the cavity under hydrostatic load is a circle of radius R i . This hydrostatic load results from a radial displacement U C prescribed on the solid at a distance R e of the center. The loading normal to the idealized grain boundary is realized by superposing to the radial displacement a vertical, negative displacement, proportional to the distance to the horizontal axis and to the positive scalar U D . This scalar is the maximum displacement magnitude prescribed and is a function of time, as shown in Fig. 13b. After a rise time 1 , the scalar U D is either kept constant or brought back to zero between times 2 and 3 . This cycle is typical of a test in the laboratory. Note that in view of the symmetry of the loading and geometry, only a quarter of the domain is studied. Material, geometrical and loading parameters are found in Table 3 if they di er from the values provided in Table 2. In this subsection the uid phase is considered as an open system so that the chemical potential is constant in time. The case of a closed system is considered in the next subsection.

Results for problem 2, a constant chemical potential and a loading held constant after time 1 , are presented in Figs. 141516. The morphological evolution of the cylindrical channel shares several features with problem 1. This evolution allows one to distinguish two stages. First, a di use ovalization of the initial circular shape sets in and, second, a rapid dissolution takes place in a narrow region crossing the idealized grain boundary (Fig. 14a). The position of the solid-uid interface is plotted on that ÿgure for 10 times deÿned in Fig. 15 which presents the time histories of the non-material velocity of the point of intersection between the interface and the horizontal axis (the idealized grain boundary). The velocity of localized dissolution is thus orders of magnitude larger than the initial evolution rate. This rapid dissolution takes place along the idealized grain boundary providing a mechanism for the uid to penetrate the dry contact region. This undercutting of the dry contact or island leads to the creation of a uid ÿlm which, at time t 9 , has a thickness 4 times the characteristic length (inset of Fig. 14a). It is thus surface tension which appears to control the ÿlm thickness in our simulation. The undercutting of the contact takes place after a latency period necessary for the build up of the stress concentration. As in problem 1, this latency period is shortened by increasing the normal load, as is seen for the solid curve in Fig. 15 obtained for a U D of 4.5. The distribution of stress in the region ahead of the dissolution front is presented in Fig. 16 in the form of 26 isocontours of equivalent Von Mises stress ranging from 1.5 to 20 (×10 MPa) and obtained at time t 9 . Note, as in problem 1, the typical stress distribution of a fracture mechanics problem in the region ahead of the undercutting. The stress concentration observed is responsible for the large Helmholtz free energy at the tip and the increasing velocity at which the undercutting of the dry contact region takes place.

=kg R = R 18 Initial diluent phase density 0 D kg=m 3 R 0.302 Diluent isothermal compressibility Ä D Pa -1 -1 R 9:869 × 10 -3 Initial solute concentration C 0 S - - 1 0 -3 Initial internal radius R i m L R 500 External radius R e m L R 5000 Constant displacement U C m L R -0:951 Maximum vertical displ. U D m L R -2:
An important issue arising with the interpretation of experimental procedure is the nature and the extent of possible changes in the microstructure during unloading and prior to the drying of a grain contact necessary for the observation. To assess the importance of the duration of unloading, we have conducted two unloadings starting at time b1 and c1 deÿned by the intersections of the dotted and dashed curves with the solid curve in Fig. 15. Unloading occurred over a period 3 -2 which is for the ÿrst and second unloading 330 times the characteristic time, respectively (see Table 3). The ÿrst unloading takes place before the localized dissolution (during stage 1) and the ovalization is reversed towards the equilibrium shape marked by a dotted curve in Fig. 14b. This reverse process is slow as can be judged by the times b 1 and b 2 , deÿned in Fig. 15, at which the two proÿles are plotted in Fig. 14b. Furthermore, the initial equilibrium shape is not stable, as discussed in [START_REF] Leroy | Geometrical evolution of stressed and curved solid-uid phase boundaries. Part 2-Stability of cylindrical pores[END_REF], and the cavity keeps on closing at a very slow rate, as can be seen by the times provided in that ÿgure. Note that the solid-uid interface of the closing tube is not circular and retains the elliptical shape inherited during loading.

The second simulation of unloading is initiated during stage 2 at time c 1 t 5 . The localized dissolution is already well deÿned (see Fig. 14a). The unloading can then also be described as a two-stage process. The ÿrst stage lasts as long as the duration of the unloading and is characterized by a healing of the newly formed ÿlm due to a rapid deposition in that region. This stage is recognized in Fig. 15 by the change in sign of the velocity of the tip of the dissolution pit to a value which is, however, of same magnitude as the dissolution velocity. Again, this change in non-material velocity occurs during the unloading 3 -2 . The second stage of the closing of the newly formed thin ÿlm occurs after time c 3 and corresponds to a rather di use morphology change over the whole interface. The rate at which the closing of the cavity occurs is then identical to the one found in the ÿrst simulation of unloading. Furthermore, by comparing the shapes of the solid-uid interface in Fig. 14c, we observe that the memory of the existence of the thin ÿlm is completely lost in the long term.

The rapid closure of the newly formed ÿlm during stage 1 of unloading deserves some comments. It can be explained by looking at the expression for the driving force in Eq. ( 12) and at the isocontours of equivalent Von Mises stress prior to unloading, which are close to those displayed in Fig. 16. The Von Mises stress is then of the order of 50 at the tip of the newly formed ÿlm and the Helmholtz free energy ( s 2 =E) is of order 1, in view of the material properties given in Table 2. The Helmholtz free energy is the leading term in the driving force and this value of 1 is consistent with the non-material velocity found at time t 5 in Fig. 15 according to the kinetic law (11). Upon unloading, the stress distribution in the region of the cusp tip changes drastically. If the ÿnal loading could be described by an external pressure identical to the internal pressure, a constant state of hydrostatic stress would prevail in the solid phase for any solid-uid interface morphology. This result is of course valid in the absence of surface tension which modiÿes the local normal stress vector acting on the solid-uid interface depending on the curvature, in accord with Laplace's boundary conditions (4). To document the e ect of curvature on the stress distribution, 26 isocontours of equivalent Von Mises stress are plotted in Fig. 17 at time c 3 . This time is just after the peak in non-material velocity during unloading (Fig. 15). The shape of the interface has not yet been modiÿed from time t 5 but the stress distribution is already very di erent from the one described above at time t 9 . Surface tension results in a positive traction at the tip of the newly formed ÿlm which reduces the value of the equivalent stress compared to the inner part of the solid phase. It is now at the tip of the ÿlm that the stress magnitude is minimum (marked by a minus sign). The maximum stress is found away from the dissolution region and is marked by a plus sign in Fig. 17. Considering now the driving force for interface migration, the Helmholtz free energy at time c 3 is 1=1000 the value taken at time t 9 . The dominant contribution to the driving force is now due to surface tension (of order 2H ). The radius of curvature being of order 1, see Figs. 16 and17, the positive value of the non-material velocity at time c 3 of order 0.3 is thus well explained by the kinetic law (11) and the driving force (12).

The cylindrical tube as a closed system

The last set of results to be presented concerns problem 2 (Fig. 13) for a closed system: the compressible uid phase is isolated within the grain boundary and the diluent mass is preserved. Di usion is fast within the uid which is assumed homogeneous at all times. The pore pressure and solute concentration are thus computed at every time step as described in Section 2 and in Appendix A. The 26 isocontours (a-z) range from 1.10 to 1.25 (×10 MPa). Note that, due to the action of surface tension, it is not the maximum, marked by a plus sign, but the minimum, found at the position of the minus sign, which is at the grain boundary contact. The L in the top right corner provides the length scale: 20 (×10 -8 m).

To understand the non-linear evolution of the pressure, the solute concentration and thus of the solute chemical potential, it is ÿrst proposed to consider the model problem 2 (Fig. 13a) in which the loading is radial at all times. The radial displacement results from the superposition of a constant U C and a variable U D having the same history as shown in Fig. 13b. The loading being radial, the problem is axisymmetric with an analytical solution for the driving force for any value of the inner radius R i . The problem solution is nevertheless time dependent and requires the numerical scheme presented by [START_REF] Ghoussoub | Solid-uid phase transformation within grain boundaries during compaction by pressure solution[END_REF]. The results are presented in Fig. 18 where the uid pressure, solute concentration, solute chemical potential and thermodynamic force are plotted as a function of time. Material, geometrical and loading parameters have the values provided in Table 3.

The evolution of the uid pressure presented in Fig. 18a is dominated by the compressibility of the uid phase, since there are small di erences between the solid and the dotted curve, the latter being obtained in the absence of a phase transformation (no PT). This latter curve is composed of straight segments connected at the end of the rise time 1 , the beginning of unloading 2 and at the end of the loading cycle 3 . There is thus little in uence on the pressure of the variation in solute concentration presented in Fig. 18b. This solute concentration decreases by 6% during the loading cycle. Note the time delays due to the kinetic law at the beginning of the cycle and at any change in loading conditions (times 1 , 2 and 3 ). Another feature seen on those ÿrst two graphs of Fig. 18 is the constant state of pressure and solute concentration reached before unloading: the system has achieved a new equilibrium which is stable to homothetic perturbation. Such a stable equilibrium shape was not present for open systems and is due to the mechanical interaction between the solid and uid phases.

The next quantity studied during the radial loading of the cylindrical cavity is the chemical potential (Fig. 18c). The solid curve represents the variation in chemical potential with respect to the equilibrium value * S . This variation is due to two contributions according to Eq. ( 9), the ÿrst resulting from the pressure dependence of the chemical potential and the second from the variation in solute concentration. Those two contributions are presented by the dotted and dotted-dashed curves in Fig. 18c, respectively. The ÿrst contribution leads to an increase in the chemical potential during loading while the second contribution is negative and is due to the decrease in solute concentration. Note that the dotted-dashed curve exhibits the delays described above when loading is initiated or modiÿed. It is the sum of these two contributions which results in the complex variation of the chemical potential, which increases ÿrst up to time 1 (in uence of p f ) and then decreases once loading is held constant (delay due to the deferred evolution of the solute concentration). This pattern repeats itself during the second part of the loading cycle. The same complexity is found for the driving force deÿned in Eq. ( 12) and plotted in Fig. 18d as a solid curve. This evolution is dominated by the variation of the uid pressure and the chemical potential. This conclusion is reached by comparing the solid, dotted and dashed curves in Fig. 18d corresponding to the driving force, the Helmholtz free energy, and the term (p f -2H )= s which is the third on the right-hand side of the deÿnition of the driving force (Eq. ( 12); note that the dimensionless -1 s is equal to 1 + tr( )). The Helmholtz free energy constitutes a small contribution to the driving force in the absence of any stress concentration. The variation of the pressure term is as important as the variation of the solute chemical potential, as it can be seen by comparing the solid curve of Fig. 18c and the dashed curve in Fig. 18d. The di erence between the latter and the former determines mostly the evolution of the driving force. Note that this driving force is close to zero during most of the interval 3 -2 , expressing the proximity of an equilibrium state. The variations of the four quantities presented in Fig. 18 are now going to be used to explain the initial morphology changes of the solid-uid interface during a loading normal to the idealized grain boundary.

The results for the evolution of the cylindrical cavity during a loading normal to the grain boundary with an isolated uid phase are presented in Figs. 1920212223. The evolution of the closed system on the time scale considered to generate the normal loading is ÿrst commented prior to the analysis of the long term evolution. The initial evolution of the pressure, the solute concentration, solute chemical potential and the maximum non-material velocity are presented in the insets of Figs. 19 and21. The application of a displacement normal to the grain boundary, on top of the radial loading constant in time, results in a compression of the uid phase and an increase in the uid pressure which is controlled by the compressibility of these two phases as for the radial loading discussed above (compare Fig. 18a and inset of Fig. 19a in which the solution in the absence of any phase transformation is indicated by a dashed curve). The evolution of the solute concentration is also similar to the one discussed previously since it decreases by 4% until the rise time 1 is reached and tends to a plateau afterwards, inset of Fig. 19b. The solute chemical potential evolves in time with the pressure and the concentration in accord with the uniform radial analysis (see inset of Fig. 19c): the chemical potential increases by 0.4 before decreasing towards a plateau once the rise time 1 is reached. The value of the driving force cannot be presented in a way similar to S since the loading results in a heterogeneous stress state around the cavity. We can nevertheless analyze the non-material velocity of the interface at the intersection of the idealized grain boundary, inset of Fig. 21. There, the velocity is positive, increasing from time zero to 1 , and decreasing afterwards as if an equilibrium state was acting as an attractor. This positive velocity signals a deposition as for a uniform radial loading. The heterogeneous stress state generated by the loading normal to the grain boundary makes itself felt mostly in the long term. However, its e ects can already be seen on the time scale of 1 from the distribution of non-material velocities over the interface presented in Fig. 20 at times 0.77, 2.30 and 2.45 times 1 . The vectors in that ÿgure are scaled with respect to the maximum magnitude in velocity which are 0.25, 1:3 × 10 -2 and 6:6 × 10 -3 at the 3 times denoted (a), (b) and (c) and deÿned in the inset of Fig. 21. The distribution in Fig. 20a shows the general deposition over the whole interface during the application of the loading normal to the grain boundary. The consequence of the stress heterogeneity on the non-material velocity is then not important. Once this loading is established, the non-material velocity decreases and the stress gradient manifests itself by a smaller non-material velocity close to the idealized grain boundary (Fig. 20b). At time 2:45 1 (Fig. 20c), this e ect is su ciently important for the velocity in the region of the grain boundary to be reversed and for the newly added material to start to dissolve. From that time onwards, the evolution of the closed system resembles more the one described in the previous section for an open system except for a major di erence which is discussed next.

The second phase of the morphology evolution occurs on a very di erent time scale which is 6 orders of magnitude larger than the rise time. This is the time scale necessary for the ovalization of the cylindrical cavity to be su cient to generate a stress concen- tration along the idealized grain boundary. This ovalization and the consequence of the stress concentration build up can be seen in Fig. 22 for 12 times which are deÿned in Fig. 21. This long term evolution is thus similar to the one discussed in the previous section: once the stress concentration is established after a latency period, there is a rapid dissolution process which provides a mechanism for the uid to wet the grain boundary and to create a uid ÿlm. The velocity of the tip of the newly formed thin ÿlm is then several orders of magnitude larger than during the latency period. The di erence between the undercutting mechanism along the idealized grain boundary for closed and open systems can now be judged by comparing the insets of Figs. 22 and14a. It can be seen that for the closed system there is a deposition along the newly formed ÿlm over the X -coordinate -650 to -750. To understand this thinning of the ÿlm, one has ÿrst to recognize that the solute concentration in the uid phase keeps on decreasing during the undercutting, (Fig. 19b). The rate of decrease correlates with the increase in tip non-material velocity. The pressure is also found to increase at a similar rate with the tip non-material velocity. This build up in pressure indicates that the volume occupied by the uid phase is decreasing during the wetting of the grain boundary. There is thus more material deposited over the interface than dissolved at the tip of the undercutting zone. The combination of the evolution of the uid pressure and the solute concentration leads to an increase in the chemical potential in time (Fig. 19c). To understand the consequences of this increase which results in a change of sign of the driving force over the solid-uid interface the various contributions to the driving force should be documented as in [START_REF] Ghoussoub | Solid-uid phase transformation within grain boundaries during compaction by pressure solution[END_REF].

The conclusion of such a detailed analysis can be comprehended by studying the distribution of equivalent Von Mises stress at 4 times labeled (d) -(g) and deÿned in Fig. 21. Fig. 23 presents 26 isocontours at each time in a quarter of the ÿgure, starting with time (d) in the upper right corner and counting anticlockwise. The isocontours range always from 1.5 to 20 (×10 MPa) for sake of comparison. It can be seen that, as the undercutting develops, the crack-like behavior results in a general unloading away from the tip. For example, the isocontour i which corresponds to a stress of ÿve intersects the cavity at a point which is shifted towards the undercutting tip as times changes from (d) to (f). At time (g), this level of stress is not found any more on the solid-uid interface within the cavity. The evolution of the stress concentration towards a crack-like stress distribution is thus responsible for the change in sign of the driving force over most of the cavity wall and, consequently, for the thinning of the newly formed ÿlm.

Conclusion

This paper has given some insight on the solid-uid phase transformation which takes place within grain boundaries of porous rocks compacting by the deformation mechanism IPS. For that purpose, two initial-and boundary-value problems have been proposed, motivated by the channel and island structure observed within grain boundaries after experiments [START_REF] Urai | Weakening of rocksalt by water during long-term creep[END_REF]Spiers et al., 1990).

The assumptions common to the two problems are that the phase transformation is controlled by a linear kinetic law which relates the mass ux to the di erence between the chemical potential of the solid phase at the interface and of the solute in the uid phase. The solid is isotropic, linear elastic and di usion occurs only within the uid phase assumed homogeneous and in chemical and mechanical equilibrium at all times. The numerical method relies on an explicit scheme and is based on the ÿnite-element method. At every time step, the state of stress is determined over the domain occupied by the solid phase and the non-material velocities of the phase transformation are computed over the solid-uid interface, thanks to a weak formulation of the kinetic law. The boundary of the solid phase is then migrated based on that velocity distribution prior to an adaptation of the ÿnite-element mesh over the solid phase. The feature of the family of ÿnite elements proposed is a C 1 continuity along the interface and a C 0 continuity in the bulk obtained by combining Hermite's and Lagrange's polynomials. This enhanced continuity on the interface was found necessary to capture properly the in uence of surface tension.

The ÿrst of the two problems is motivated by the possible existence of uid ÿlms within grain boundaries. The solid phase has initially the shape of a rectangular block and is in contact on one of its side with the uid phase of the ÿlm. The stability of the at interface is studied ÿrst by a perturbation analysis which reveals that the di erence between the interface parallel compressive stress and the uid pressure is destabilizing. Surface tension, on the contrary, is stabilizing. The combination of the two opposing factors results in the selection of a dominant perturbation with a speciÿc wavelength having the fastest rate of growth. There are two interesting outcomes from this linear stability analysis conducted on a domain of ÿnite dimensions. The ÿrst is that the limit for an inÿnitely long wavelength does not correspond to the stability predictions obtained for a homothetic mode. This is due to the presence of a singularity in the stability boundary conditions. This result could explain similar singular limits found in bifurcation and stability analyses such as the one conducted by [START_REF] Leroy | Spatial patterns and size e ects in shear zones: a hyperelastic model with higher-order gradients[END_REF] to detect patterns in shear zones. The second advantage of this linear stability analysis is the validation of the numerical scheme by comparing analytical predictions and numerical results during the initial development of the perturbations. The long term evolution of the trigonometric perturbation obtained by numerical means is described as a two-stage process. During stage 1, the morphology evolution is di use over the whole interface: the amplitude of the trigonometric perturbation increases with time if the interface is unstable. Stage 2 is characterized by the formation of a cusp or dissolution pit in the region where the initial perturbation resulted ÿrst in an increase in the solid Helmholtz free energy. The curvature at the tip of the cusp is related to the characteristic length due to surface tension at the end of the simulation. However, the rate of dissolution keeps on increasing at all times revealing the crack-like behavior at the tip of the cusp. The latency period, deÿned by the duration of stage 1, is observed to decrease with an increase in the magnitude of the initial compressive stress and in the initial amplitude of the perturbation. It was also shown that a supersaturated uid in contact with the solid does not prevent such localized dissolution to occur if the stress concentration is su cient.

The geometry of the second problem is a cross section through a typical channel observed within the grain boundary. The equilibrium cross section is circular since the excess energy associated to the solid-solid interface is disregarded. The compression which occurs typically during the experiments is modeled by prescribing a loading normal to the idealized grain boundary. This loading is ÿrst studied for an open system typical of a uid phase in contact with the remote pore uid. In that instance, the solute chemical potential and the uid pressure are assumed constant in time. The evolution of the circular cavity occurs, as in problem 1, in two phases. There is ÿrst an ovalization of the circular cross section at a very slow rate. Once the stress concentration is su cient along the idealized grain boundary, a rapid dissolution occurs in that region. This dissolution provides a mechanism for the uid to penetrate the whole grain boundary by undercutting the dry contact regions and thus to create a new uid ÿlm. The thickness of this ÿlm is governed by the curvature at the tip of the undercutting region which is of the same order of magnitude as the characteristic length of the problem due to surface tension. If the loading cycle is terminated by bringing back to zero the normal load, the evolution of the cavity is then also described in two stages. In the ÿrst stage, the newly formed ÿlm closes at a rapid rate of same magnitude as the undercutting rate. This fast rate of healing is due to surface tension which enters the driving force for phase transformation. Then, in the second stage, the cavity retains an elliptical shape and closes down in a homothetic way at a very slow rate compared to the ÿrst stage. The cavity analysis is repeated for a closed system characteristic of a uid phase isolated within the grain boundary. The solute concentration and the uid pressure are then two of the problem unknowns. The evolution of the cylindrical cavity as a closed system is similar to the one described for an open system except for a major di erence: the stress concentration at the tip of the undercutting region results in a general decrease in the stress magnitude around the cavity. This variation in stress leads to a change in sign of the driving force away from the tip. Consequently, there is deposition over the whole domain except in the tip region, the solute concentration drops, the volume of the cavity decreases and the uid pressure increases. Moreover, the thickness of the newly formed ÿlm decreases in time.

The relevance of these numerical results for a better understanding of the action of IPS within a grain boundary should now be discussed. It should be stressed that the wavy perturbations of problem 1 are observed experimentally on K-alum salt in the laboratory [START_REF] Den Brok | The e ect of elastic strain on the microstructure of free surfaces of stressed minerals in contact with an aqueous solution[END_REF]. Furthermore, the cusps observed in problem 1 are reminiscent of the observations made within the grain boundary of rocks which deformed by IPS in nature [START_REF] Welton | SEM Petrology Atlas. Methods in Exploration Series[END_REF]. The long latency period prior to the localized dissolution, obtained with a data set typical of quartz grains at 200

• C, deserves some comments. During that period the grain boundary is partly dry and thus IPS inactive. This long period is in line with the di culty to produce IPS with quartz in the laboratory. A better quantitative assessment of the latency period prior to the action of IPS would require a characterization of the grain boundary roughness to estimate the relevant stress concentration factor. Once this stress concentration is su cient, we observe a fast undercutting of the dry regions. It seems that our numerical results are the ÿrst proof that a uid can wet a dry contact region by marginal dissolution resulting in the formation of a new uid ÿlm. If loading is interrupted, the newly formed ÿlms close partly during the lapse of the unloading. Consequently, part of the structure operating during IPS is likely to be lost when the physical observation is conducted after the experiment. There are, however, no contradictions between the observation of a channel and island structure after the experiment and the presence of thin ÿlms during the same experiment. Furthermore, extrapolating the stability analysis of our ÿrst problem to the thin ÿlms created in problem 2, it appears that this thin ÿlm is not stable if the interface parallel stress di ers from the pore pressure. One should thus expect the thin ÿlm morphology to evolve toward a corrugated surface which could end up with the creation of new solid to solid contacts. This thinning of the ÿlm is found here to be enhanced in closed systems. This possible continuous reorganization in time of the grain boundary points to the direction of a dynamic internal structure within the grain boundary which has been often invoked in the literature.

Eq. (B.1) following the procedure used for Lagrange elements in the literature (e.g., [START_REF] Hughes | The Finite Element Method[END_REF].

The family of isoparametric and compatible elements described above is complete. A uniform gradient over the element is inferred by nodal values calculated from any constant gradient ÿeld [START_REF] Ghoussoub | Solid-uid phase transformation within grain boundaries during compaction by pressure solution[END_REF]. These conditions are su cient to ensure the convergence of this class of ÿnite elements for elliptic problems including the mechanical equilibrium (13) and the mesh adaptation problems. normal displacement and zero shear stress conditions reveal that the constants a are real and that a 1 = a 2 = 0. Along the solid-uid interface, the linearization of Laplace equations ( 4), together with n = -i e x ; t = i e y and 2 Ĥ = -2 (according to simple surface geometry results; [START_REF] Aris | Vectors, Tensors and the Basic Equations of Fluid Mechanics[END_REF], provides two conditions for a 3 and a 0 which, once inverted, lead to a 0 a 3 = S + HC HS -S -C

1 p f + 0 × 1 CS + H : (C.4)
The letters C and S in Eq. (C.4) stand for cosh( H ) and sinh( H ), respectively. The last boundary condition to be inspected is along A B and depends on whether force or displacement is prescribed on that surface. If the normal displacement is prescribed, its perturbation is zero and the wavenumber has to be equal to k =L ( with k ∈ N) while the no-shear stress conditions is then satisÿed by the previous ÿndings on the constants. For force control, the total force on the boundary remains constant constraining the perturbation in stress by This condition leads to the same equation as the second equation obtained from the linearized Laplace conditions. The second condition for force control is a no-shear stress conditions which is found to be satisÿed for the same reason mentioned above if = L is set to any positive, non-zero integer k.

Further calculations based on the linearization of the kinetic law (11), which are not discussed here for sake of brevity, leads to the following expression for the stability exponent:

= -2 1 - ( 0 + p f ) 2G G 1 ( H ) + ( 0 + p f ) 2 (1 -) G C 2 CS + H with G 1 ( H ) = H -(1 -2 )CS CS + H ; (C.6)
which is found to be real. The function G 1 of H (a multiple of H =L) takes values between 1=2 and 1. Consequently, the factor multiplying 2 in Eq. (C.6) should be set to ÿrst order to 1 to be consistent with the small-strain hypothesis. This simpliÿcation leads to expression (24) in the main text for the stability exponent. This expression (24) for the stability exponent is commented in the main text and we only focus in the rest of this appendix on the limit taken by as tends to zero. This limit has to be considered for inÿnitely long rectangular block since remains equal to k =L to satisfy the boundary conditions postulated for the linear stability analysis. Taking this limit in Eq. (C.6), one ÿnds

lim →0 = ( 0 + p f ) 2 (1 -) G 1 2H : (C.7)
This limit is positive or zero and is independent of the type of boundary condition adopted along A B . This is of course in contradiction with the analysis of the fundamental solution commented with the support of Fig. 5. There, it is seen that the fundamental solution is neutrally stable for displacement control and always unstable for force control. Furthermore, a rapid calculation on the fundamental solution for force control [START_REF] Ghoussoub | Solid-uid phase transformation within grain boundaries during compaction by pressure solution[END_REF] reveals that the stability exponent for a homothetic perturbation, a perturbation in H which preserves the rectangular shape of the block, is identical to the limit found in Eq. (C.7). This discrepancy for displacement control is clariÿed by looking again at the stability analysis for wavy perturbations noting ÿrst that the stress perturbation should be zero for a displacement control test and a perturbation of inÿnitely long wavelength. The constants a 0 and a 3 should then be at least of order 1 in . However, it turns out that the system of equations for (a 0 ; a 3 ) obtained from Laplace conditions is singular for a zero . Its solution for a 3 presented in Eq. (C.4) has a leading term in ( H ) -1 which governs in the limit of zero H . It is this singular behavior which results in the di erence between the stability exponent in the long-wavelength limit and for homothetic change in thickness.

The structure of the stability analysis just presented is common to many bifurcation and stability problems. We have seen here that a singular limit could explain a discrepancy between homothetic and long-wavelength mode. It is worth noting that a similar discrepancy was found by [START_REF] Leroy | Spatial patterns and size e ects in shear zones: a hyperelastic model with higher-order gradients[END_REF] for the same type of limit in a 2D shear band analysis. The present simple calculation could provide a ÿrst insight on the origin of their results.

Fig. 1 .

 1 Fig. 1. The ÿnite-element discretization. The C 1 discretization of the interface (a) is realized with Hermite polynomials. This interpolation is coupled to Lagrange polynomials to construct a ÿnite-element with C 0 continuity in the bulk and C 1 on the interface (b). The local node numbering of the element is presented in the computational space in (c).

Fig. 2 .

 2 Fig. 2. Transmission optical micrograph showing uid inclusions within grain boundaries in experimentally deformed rocksalt. The salt is synthetically produced and contains about 50 ppm of water. It was deformed at a temperature of 200 • C, a conÿning pressure of 50 MPa and a constant strain rate of 5 × 10 -7 s -1 . The tubular inclusions are typically 3 m in diameter. Results were obtained in the HPT Laboratory, Faculty of Sciences, Utrecht University, The Netherlands and presented with the kind permission of Prof. C.J. Spiers.

Fig. 3 .

 3 Fig. 3. Schematic diagram of a grain boundary and of two di erent features of its internal structure. They lead to the proposition of two initial-and boundary-value problems. The ÿrst problem (1) consists of a rectangular block sustaining a plane-strain compression and exposed on its top side to the uid phase. The second problem geometry (2) is a cylindrical cavity, a quarter of the structure being studied for reason of symmetry.

Fig. 4 .

 4 Fig. 4. The ÿrst model problem consists of a rectangular block in contact on its top side with the uid phase. The wall A B is displaced uniformly to produce a compression by either keeping the force F constant or prescribing the displacement U . The boundary conditions on the two sides AB and AA consist of no normal displacement and, in the absence of friction, of a zero shear stress. Laplace conditions (Eq. (4) in the text) are imposed on the solid-uid interface.

Fig. 5 .

 5 Fig.5. Phase diagram for the fundamental solution of problem 1 deÿned by the rectangular block of thickness H sketched in the previous ÿgure. The rate of dissolution or precipitation is independent of H in a displacement control test (DC), as shown by the dotted lines for the special case of 0 = -p f . For the force control case (solid curves), an equilibrium (solid dot) to exist requires the dimensionless chemical potential to be larger than the pore pressure. Note that this equilibrium is unstable for a homothetic perturbation.

Fig. 7 .

 7 Fig. 7. The evolution in time of the perturbation size is analyzed by plotting the d 1 distance function (surface between the current solid-uid interface and the initial at perfect interface). This measure is normalized by the value of d 1 at time zero. The perturbation is a long-wavelength mode with = =L and L = H . The agreement between the numerical solutions (solid curves) and the analytical predictions based on the linear stability analysis (dotted lines) validates the numerical scheme proposed.

Fig. 11 .

 11 Fig. 11. The mesh and isocontours of equivalent Von Mises stress are presented in the tip region of the dissolution pit. The 26 isocontours (a-z) range from 15 to 50. The stress is approximately 150 at the tip (×10 MPa). The L in the center provides the scale in the two directions (ratio 1 : 2) and the horizontal bar has for length 1.8 (×10 -8 m).

Fig. 13 .

 13 Fig. 13. Model problem 2 consists of a quarter of a cylindrical cavity with the geometry and loading shown in (a). The loading on the external radius R e at equilibrium consists of a radial displacement U D which remains unchanged with time. To produce a loading normal to the idealized grain boundary, a vertical displacement proportional to the vertical coordinate and having for maximum U D is applied. The maximum value U D has the history shown in (b) and is either kept constant after a rise time 1 or brought back to zero between 2 and 3 .

Fig. 14 .

 14 Fig. 14. Morphology evolution of the cavity, which equilibrium shape is circular, during a loading cycle normal to the idealized grain boundary (horizontal axis). In (a), the loading is held constant after a short rise time and dissolution is localized along the grain boundary. In (b) and (c), the loading is brought back to the equilibrium condition in a short time lapse. The various times are identiÿed in Fig.15.

Fig. 15 .

 15 Fig. 15. Evolution with time of the non-material velocity at the point on the idealized grain boundary. The solid lines are time histories for a compressive loading normal to the grain boundary. The dotted and dashed lines are two velocity histories as the loading is brought back to the initial hydrostatic loading condition. A time of 1000 and a velocity of 1 represent 20.3 yr and 0.5 m=yr, respectively.

Fig. 16 .

 16 Fig. 16. The mesh and 26 isocontours of equivalent Von Mises stress in the region ahead of the localized dissolution at time t 9 deÿned in Fig. 15. The isocontours (a-z) range from 1.5 to 20 (×10 MPa). The L on the right provides the scales in the two directions (ratio 2 : 1) with the horizontal bar having for length 33.3 (×10 -8 m).

Fig. 17 .

 17 Fig.17. Isocontours of equivalent Von Mises stress are plotted in the zone close to the grain boundary where localized dissolution took place before unloading. The 26 isocontours (a-z) range from 1.10 to 1.25 (×10 MPa). Note that, due to the action of surface tension, it is not the maximum, marked by a plus sign, but the minimum, found at the position of the minus sign, which is at the grain boundary contact. The L in the top right corner provides the length scale: 20 (×10 -8 m).

Fig. 18 .

 18 Fig. 18. Evolution with time of the uid pressure (a), the solute concentration (b), the solute chemical potential (c) and the thermodynamic force (d), during a radial loading cycle on the cylindrical cavity of Fig. 13. The radial displacement is the constant U C plus the scalar U D , the history of which is shown in Fig. 13b.

Fig. 19 .

 19 Fig. 19. Evolution with time of the uid pressure (a), the solute concentration (b), and the variation of the solute chemical potential (c) during the loading of the cylindrical cavity shown in Fig. 13. The uid phase is assumed isolated. In insets, the evolution of the three quantities of interest is shown on the time scale of the rise time.

Fig. 20 .

 20 Fig.20. The non-material velocity on the solid-uid interface at three di erent times, 0:77; 2:30 and 2.45 times the rise time 1 in (a), (b) and (c), respectively. The size of the vectors are normalized by the maximum velocity magnitude which is 0:25; 1:3 × 10 -2 and 6:6 × 10 -3 at the three times.

Fig. 21 .

 21 Fig. 21. History of the non-material velocity at the point of the solid-uid interface in contact with the idealized grain boundary. In the inset, this history is presented on the time scale of the rise time.

Fig. 22 .

 22 Fig. 22. The position of the interface for the problem of a cylindrical cavity of initially circular cross section and for a closed system, at 12 deÿned on the previous ÿgure.

Fig. 23 .

 23 Fig. 23. 26 isocontours of equivalent Von Mises stress ranging from 1.5 to 12 (×10 MPa) at four times labeled (d) -(g) and deÿned in Fig. 21.

H 0 ˆ

 0 xx (y) dy + ( 0 + p f ) = 0: (C.5)

Table 2

 2 Physical constants, geometrical parameters and loading conditions used for the analysis of the rectangular block, model problem 1

	Parameter	Unit	Reference	Dimensionless
			quantity	value
	Young's modulus E	Pa	R	8200
	Elasticity modulus in shear G	Pa	R	3000
	Surface tension	Pa m	R	1
	Constant solute chemical potential S	Pa m 3 =kg	R = R	[3.012; 3.04]
	Solid phase density 0	kg=m 3	R	1
	Kinetics constant L	kg s=m 4	L R	1
	Initial thickness H = L	m	L R	100, 500 or 1000
	Constant or initial uid pressure, p f or p 0 f Initial stress 0	Pa Pa	R	3

R

[ -22; -18]

Table 3

 3 Physical constants, geometrical parameters and loading conditions used for the analysis of model problem 2, if di erent from those presented for problem 1 in Table2

	Parameter	Unit	Reference	Dimensionless
			quantity	value
	Constant solute chemical potential S	Pa m 3 =kg	R = R	2.997
	Constant scalar RT=M S	Pa m 3

I (t)

Acknowledgements

Many of the ideas presented here are the fruits of numerous discussions with W.K. Heidug (Shell Research, RTS, The Netherlands), F.K. Lehner (University of Salzburg, Austria) and the researchers of the HPT Laboratory (Utrecht University, The Netherlands). The ÿrst developments of the numerical algorithm for pressure solution were realized with L. Bodin (University of Metz, France) and P. Onck (T.U. Delft, The Netherlands) during their visit to Shell Research. J. Ghoussoub has beneÿted from a fellowship of Ecole Nationale des Ponts et Chaussà ees during her doctorate studies and the project was partly ÿnanced by a grant from the French Group of Research (Geomechanics of Deep Rocks) from 1995 to 1998. This support is gratefully acknowledged.

Appendix A. The constraint for a closed system

The objective of this appendix is to derive the constraint between the volume occupied by the uid phase, the solute mass and the uid pressure in a closed system with a ÿxed diluent mass.

The solution is dilute and the solute assumed to behave as an ideal gas sustaining the pressure p f . Consequently, the volume occupied by the solute is

from a reference state characterized by the pressure p 0 f , the volume V 0 S and the mass m 0 S . The diluent, say water, is assumed to have a constant isothermal compressibility Ä D such that the volume it occupies reads

for a reference volume V 0 D . The total volume occupied by the uid phase is the sum of the solute (A.1) and the diluent (A.2) volume as for an ideal mixture and is presented in Eq. ( 10) in the main text. It is a function of the pore uid pressure and solute mass.

The evolution of the solute chemical potential in a system constrained by a constant mass of diluent requires the knowledge of the pore pressure and the solute concentration. The current solute mass and the cavity volume are found by integrating in time ṁS = -

prior to the computation of the solute density ( S = m S =V S ) and concentration (C S = m S =(m S + m 0 D )). The uid pressure is then found by solving Eq. ( 10) and the solute chemical potential from Eq. ( 9).

Appendix B. A class of ÿnite elements for IPS

In this appendix one deÿnes the shape functions for the class of elements employed in the numerical simulations. The convergence property of the elements is also brie y analyzed.

The shape functions for the basic four-noded element (plus two tangents) having a C 1 continuity along the side with 2 equal to +1 result from the following products of Lagrange and Hermite polynomials:

where the local numbering is deÿned in Fig. 1. Additional nodes (7-10) can be included with a construction of their shape functions and modiÿcations of the ones presented in

Appendix C. Linear stability of the rectangular block

This appendix provides a summary of the derivation of the stability exponent , found in the main text as Eq. ( 24), which characterizes the initial time evolution of trigonometric perturbations added to the otherwise at solid-uid interface of the rectangular block presented in Fig. 4. Details of the analysis are found in [START_REF] Ghoussoub | Solid-uid phase transformation within grain boundaries during compaction by pressure solution[END_REF]. This linear analysis di ers from the ones published so far by the account of the ÿnite size of the domain. Results in the long-wavelength limit are di erent from those obtained for homothetic growth (the geometry then remains rectangular) due to a singularity which is discussed. The analysis of this singular limit could shed light on other stability and bifurcation problems having similar limits [START_REF] Leroy | Spatial patterns and size e ects in shear zones: a hyperelastic model with higher-order gradients[END_REF] in view of the common structure of the analyses.

The starting point of the linear stability analysis is the proposition that the interface position is modiÿed by an inÿnitesimal amount according to the real part of

in which H , i and are the thickness of the rectangular block at equilibrium, the pure imaginary number and the wavenumber of the perturbation, respectively. If the solution of the linearized boundary-value problem proves that the real part of is negative for all admissible , then the solid-uid interface is stable. A positive real part signals the onset of instability at the wavenumber .

The inÿnitesimal change in the position of the interface results in any ÿeld variable A(x; y) to be perturbed by a function A(x; y) which is decomposed in normal mode as in Eq. (C.1). The normalized function A(x; y) takes the form

in terms of the unknown  (y), a function of the y-coordinate only.

The ÿeld equilibrium equations over the block of perturbed geometry are satisÿed by the introduction of the Airy stress function (x; y). Compatibility equations are then found to require the Laplacian of (x; y) to equal zero which dictates that the function ˆ (y) has the following structure: ˆ (y) = (a 0 + a 1 y) cosh( y) + (a 2 + a 3 y) sinh( y);

(C.3) in which a ( = 0; : : : ; 3) are four complex constants. These four constants are determined by the boundary conditions over the 4 sides of the rectangular block shown in Fig. 4. Along the sides AA and AB, the zero