Microscale THMC Modeling of Pressure Solution in Salt Rock: Impacts of Geometry and Temperature

Pressure solution, a mechanism that involves tight coupling between the geometry and thermal-hydro-mechanical-chemical (THMC) processes, plays an important role in diagenesis. In this study, we make the first attempt to conduct microscale THMC modeling to understand and quantify the impacts of geometry and temperature on pressure solution, taking natural salt rock as an example. This modeling capability is achieved by expanding a novel MC code that we developed previously (Hu et al. J Geophys Res: Solid Earth 126:e2021JB023112, 2021) to include temperature effects. We first conduct a simulation of an example that involves a single brine inclusion within a single halite grain and find that the temperature impact is limited for that case. We then extract geometry from an image of a natural salt rock and conduct simulations with different cases: (A) only temperature and no stress, (B) only stress and no temperature, and (C) with both stress and temperature. These different cases result in quite different phenomena. In case A, dissolution and precipitation occur across the entire system due to isolated pore space reaching a localized mass balance between dissolution, precipitation, and diffusion. In case B, intense geometric features (e.g., major asperities, inclusions) in one area undergo stress concentration, thus dominating pressure solution in that area. In case C, pressure solution is spread out at contacting highly stressed geometric features close to the hotter side. We conclude that geometric features dominate stress distribution, thus dominating pressure solution in a natural salt rock that may be affected by the temperature if a sufficient temperature gradient is applied. First microscale THMC model for analyzing pressure solution that considers realistic geometry and complex multiphysics. Geometric features dominate pressure solution in a natural sedimentary rock because they dominate the stress distribution. Pressure solution increases the connectivity in a system, thus potentially increasing dissolution over a large area where stress is sufficiently high. Temperature gradients may have a limited impact on pressure solution, but can change the system behavior by changing the magnitudes of pressure solution in localized areas. First microscale THMC model for analyzing pressure solution that considers realistic geometry and complex multiphysics. Geometric features dominate pressure solution in a natural sedimentary rock because they dominate the stress distribution. Pressure solution increases the connectivity in a system, thus potentially increasing dissolution over a large area where stress is sufficiently high. Temperature gradients may have a limited impact on pressure solution, but can change the system behavior by changing the magnitudes of pressure solution in localized areas.

• First microscale THMC model for analyzing pressure solution that considers realistic geometry and complex multiphysics.• Geometric features dominate pressure solution in a natural sedimentary rock because they dominate the stress distribution.• Pressure solution increases the connectivity in a system, thus potentially increasing dissolution over a large area where stress is sufficiently high.• Temperature gradients may have a limited impact on pressure solution, but can change the system behavior by changing the magnitudes of pressure solution in localized areas.

Introduction
Pressure solution plays an important role in the Earth's sedimentary rocks (Sorby 1863; Rutter 1976Rutter , 1983;;Tada and Siever 1989;Spiers et al. 1990;Gratier et al. 2013).Pressure solution, together with other granular processes such as dislocation and micro fracturing, affects the long-term time-dependent deformation of the Earth systems (Urai et al. 1986;Olivella and Gens 2002;Hu et al. 2021).These long-lived evolving Earth systems are fluid-filled, heated, and stressed, and they consist of different types of minerals that may be dissolved or grow as a result of coupled thermalhydro-mechanical-chemical (THMC) processes.A salt rock is an example of sedimentary rock.Pressure solution takes place at the grain scale, and affects the creep of salt at larger scales (Urai et al. 1986;Hu et al. 2021).The creep of salt at larger scales is a key control on the long-term safety of engineering activities such as nuclear waste disposal and oil/ gas storage in salt rocks (Urai et al. 1986;Kuhlman et al., 2020).In turn, the mechanical, chemical, and thermal loading/unloading that are applied by humans in those engineering activities such as excavation, storing chemically active fluid (such as hydrogen), and disposing heat-generating waste packages may affect the pressure solution, potentially changing the creep behavior of the natural salt rocks.Thus, fully understanding the role of coupled THMC processes that contribute to pressure solution in natural settings is essential for understanding the formation and evolution of the sedimentary rocks in our planet Earth's evolution, and for effectively making use of these systems for safe engineering disposal and storage.A large number of experimental studies have been conducted to understand the microscale processes including pressure solution for different rock materials that include rock salts (Urai et al. 1986;Spiers et al. 1988Spiers et al. , 1990;;Peach et al. 2001;Ter Heege et al. 2005;Urai and Spiers 2007;Zhang et al. 2007;Linckens et al. 2016;Mills et al. 2018;Kuhlman et al. 2020), gypsum (De Meer et al. 2000), calcite (Zhang and Spiers 2005), and quartz (van Noort et al. 2008).These studies have greatly advanced our understanding of pressure solution and creep from micro to core scales.In addition, salt samples that were obtained from different fields have been analyzed in the laboratory tests with micrographs for using microstructural signatures to explain deformation and recrystallization mechanisms (Schléder and Urai 2007;Desbois et al. 2010;Schoenherr et al. 2010).Despite the importance of these experimental and characterization studies, it remains critical to have an effective quantitative toolset that can be used to quantify each of the mechanisms that contribute to different components of THMC coupling.This is essential for predicting the long-term behavior that cannot be investigated by laboratory tests alone.
With the development of theoretical and numerical toolsets, a number of models have been developed for quantifying pressure solution (Renard et al. 1997;Gundersen et al. 2002;Yasuhara et al. 2003;Yasuhara and Elsworth 2004;Pluymakers and Spiers 2015;van den Ende et al. 2018van den Ende et al. , 2019)).Due to the grand challenges associated with contact dynamics among deformable Earth material bodies which may have arbitrary geometry and structures, it was assumed in those models that material bodies are either spherical or regular.Or the changes in where and how contacts occur between the material bodies are not captured, for example, the transition from being not in contact to being in contact, or the transition from being in contact to sliding along the surfaces.It is further complicated by the fact that pressure solution involves diffusion and chemical reaction, thus requiring a rigorous treatment of coupling of deformation, contact dynamics, flow transport, chemical reaction, and evolving geometry.To address these challenges, Hu et al. (2021) developed a first-of-its-kind mechanical-chemical (MC) model that rigorously accounted for the microscale MC coupling with realistic geometry and multiphysics considering grain relocation, microfracturing, and pressure solution.This MC is realized by linking a mechanical code named Numerical Manifold Method (NMM, Hu and Rutqvist 2020b) to the reactive transport code named CrunchFlow (Steefel et al. 2015a).
The numerical manifold method (NMM, Shi 1992Shi , 1996) is a method for analyzing the mechanics of both continuous and discontinuous media.In NMM, two different types of meshes are defined: i.e., mathematical and physical meshes.The mathematical meshes cover a material domain of interest, and the physical meshes are divided from the mathematical meshes by boundaries and discontinuities.Independent functions are defined in different physical meshes and the weighted average of these physical mesh functions defines the approximation of an element-the overlapping area of physical meshes.Based on this, both continuous and discontinuous processes can be solved by flexibly defining physical mesh functions.To date, NMM has been widely applied for analyzing crack growth (Ning et al. 2011;Zheng and Xu 2014;Zheng et al. 2020), hydraulic fracturing (Yang et al. 2018), microscale to macroscale rock failure (Wu et al. 2017(Wu et al. , 2018)), two-phase flow in fractured media (Ma et al. 2017), coupled flow and mechanics processes in fractured media (Wu et al. 2019;Sun et al. 2022), and grouting in fractured tunnels (Xu et al. 2021).Benefiting from these groundbreaking fundamentals, the authors have developed a number of models for analyzing flow and fully coupled hydro-mechanical processes of fractured, granular and porous media at different scales (Hu et al., 2015a(Hu et al., , 2015b;;Hu et al. 2016Hu et al. , 2017aHu et al. , 2017b;;2020a, 2020b;Wang et al. 2016).
The CrunchFlow family of codes is Open Source software for calculating chemical reaction and transport ranging from nanoscale to pore scale (Steefel et al. 2015a).With the advanced capabilities of handling reactive transport at the microscale, CrunchFlow explicitly solves for the mass balance of each chemical species affected by transport via advection and diffusion, and by chemical reactions that are described with a variety of rate laws based on thermodynamics and/or kinetics.
In this study, we expand the capabilities of the MC model to consider the temperature effects on pressure solution.We will first introduce a new definition of pressure solution in Sect. 2 with an emphasis on geometric effects and its interplay with THMC coupling.Then we will briefly introduce our modeling capabilities including the numerical formulation, coupling physics, and coupling geometry in Sect.3. In Sect.4, we will present a few scenarios of numerical simulations.First, we will fully investigate temperature effects on a single halite grain with a brine inclusion.Then we will conduct analyses and comparisons of modeling results of pressure solution in a natural salt rock by with or without considering stress or temperature to provide a quantitative description of how much geometry and temperature affect pressure solution in natural salt rock.

Pressure Solution: THMC Coupling with Geometry of Geomaterials
It has been well-established that pressure solution involves the processes of (1) dissolution of solid mass at the boundaries of grains which undergo high stress (exceeding a certain threshold), (2) diffusive transport of the dissolved solid mass in the pore space, and (3) precipitation of the solid mass when the concentration reaches or exceeds the solubility which is affected by the stress distribution (Rutter 1976;Spiers et al. 1990;Renard et al. 1997;Renard et al. 2001;Yasuhara et al. 2003;Yasuhara and Elsworth 2004;Hu et al. 2021).It is also important to note that the solute concentration is affected by the diffusive transport and chemical reaction, while the solubility is affected by stress, pore pressure, and temperature.However, there has been insufficient attention to the fact that pressure solution involves a tight coupling of THMC processes with the geometry of a granular system-because the stress and contact dynamics are governed by the geometry of the natural grains, which may have arbitrary shapes and structures.
Figure 1 shows the schematic and geometry of pressure solution at two rough solid interfaces.These two interfaces can be portions of two boundaries of two grains at a scale where the roughness of each grain boundary is visible.These two interfaces can be two sides of a rough fracture with non-evenly distributed asperities.The purpose of this figure is to be general enough to include all the possible geometry of two interfaces of geo-materials.Furthermore, we classify geometric features into the following categories: first-order geometric features that include the outlines of a grain that depicts its shape and major corners/asperities, second-order geometric features that include smaller asperities that describe the roughness of a boundary, and Nth-order geometric features that become visible (or detectable) if we zoom into a smaller scale of each grain.We make the following hypothesis: Depending on the stages of diagenesis (i.e., from grain aggregates to sedimentary rocks), first-order, second-order, and Nth-order geometric features sequentially play key roles in dominating contact dynamics, contact stress, and pressure solution in a sedimentary system.
This hypothesis can be partly demonstrated by our previous work (Hu et al. 2021) where we investigated the compaction of a loosely packed halite aggregate.We found that the sharp corners (i.e., defined as first-order geometric features) control contact dynamics (involving displacements and deformation), contact stress distribution, and pressure solution.During the processes of compaction and deformation, Fig. 1 The schematic and geometry for pressure solution at two rough interfaces some sharp corners are gradually removed via pressure solution, and the solid mass is redistributed via diffusion and precipitation.In this study, we will focus on a salt rock at a later stage of diagenesis than a halite aggregate discussed in Hu et al. 2021 so as to investigate the impacts of secondorder geometric features on pressure solution.
In addition, natural salt rocks often include a number of brine inclusions.These brine inclusions are not only geometric features in mechanics that introduce stress concentration, they also provide surface areas for dissolution and precipitation to occur within a grain.And the dissolution and precipitation in a brine inclusion can be accelerated by temperature and/or stress.One example is the migration of brine inclusions as a result of a temperature gradient (Olander 1984;Caporuscio et al. 2013;Hu and Rutqvist 2020c).
Building on the previous research, in this study we aim to expand our modeling capabilities, test the aforementioned hypothesis, and answer the following scientific question: • How do the geometry and temperature quantitatively affect the pressure solution in a natural sedimentary rock such as a salt rock?
Based on NMM, we have developed multi-scale hydromechanical models and a series of computer codes subsequently for fractures and porous media with rigorous geometrical and physical representations (Hu andRutqvist 2020a, b, 2022).At micro-discontinuum scales, we developed an algorithm that (1) rigorously detects where contact occurs for geomaterial bodies with arbitrary shapes, (2) enforces contact constraints for different contact states, and (3) iterates to reach the convergence of contact states.Thus, our microscale NMM model can simulate deformation and failure of the solid grains, and dynamic changes in contact locations and contact states, which includes intergranular frictional sliding (Hu and Rutqvist 2020b).
CrunchFlow (Steefel et al. 2015a) explicitly solves for the mass balance of each chemical species affected by transport via advection and diffusion, and by chemical reactions that are described with a variety of rate laws based on thermodynamics and/or kinetics.
NMM and CrunchFlow are coupled sequentially, as shown in Fig. 2. NMM calculates the displacement, deformation, and stress and updates the geometry each time.The updated geometry as a result of displacement, deformation and stress that are calculated by NMM are then transferred to CrunchFlow.Correspondingly, CrunchFlow updates the geometry and mesh, and uses the stress components to calculate the halite solubility based on the local thermodynamic pressure (represented by the minimum principal stress) and reaction rate.In turn, chemical dissolution or precipitation calculated by CrunchFlow may lead to changes in the geometry and/or mechanical properties, and those are calculated in NMM.Because mechanical processes (e.g., deformation and movement) occur much faster than chemical reactions and fluid transport, all the components are flexibly treated as coupled or decoupled depending on the rates of processes and the time step chosen.
As mentioned in Sect.2, pressure solution often involves a tight coupling of THMC processes with the geometry of a granular system.Thus, calculating pressure solution needs to rigorously address two aspects: (1) The physics of pressure solution-i.e., how stress impacts dissolution and precipitation.We use the term "coupling physics" to refer to the mathematical and numerical formulations that we have implemented to account for the physics of pressure solution.
(2) The geometric changes during pressure solution as a result of mechanical (motion and grain deformation), and mechanical-chemical (pressure solution and precipitation) processes.We use the term "coupling geometry" to refer to the mathematical (geometric) and numerical (computational geometric) formulations that we have implemented to account for the geometry and geometric changes during pressure solution.

Coupling Physics: Stress-Dependent Solubility Formulation for CrunchFlow
Because the solubility of halite depends on both the solid and aqueous phases, we consider the thermodynamic properties for each at the same point in space within the domain.In this respect, our approach contrasts with that taken in many studies of pressure solution where dissolution and precipitation are assumed to take place at different places within the numerical sub-grid volume.Thus, we take a pore-scale approach in which local dissolution and precipitation reactions involve the solid and solution with the same thermodynamic properties (temperature and pressure/stress) at the same point in space.Precipitation, if it occurs, occurs at a different point in space within the domain then the pressure dissolution as a result of variable chemical potential conditions related to the pressure (stress) and temperature.
In this work, we use the minimum principal stress (i.e., the maximum compressive stress) to calculate the pressuredependent solubility of the halite.
The stress effect on mineral solubility can be calculated using thermodynamic pressure as the master variable.From basic thermodynamic principles, the solubility as a function of pressure is given by (Anderson 2005;Appelo et al. 2014): where ln K P,T is the equilibrium constant at the pressure of interest, ln K P=1,T is the equilibrium constant at 1 atmosphere pressure, ΔV r is the volume change of reaction cm 3 mol −1 , T is the absolute temperature ( • K) , and R is the gas constant (82.0574587 cm 3 ⋅ atm ⋅ K −1 mol −1 ).The volume change of reaction is given by the sum of the apparent partial molar volumes of the ions Na + and Cl − in solution minus the molar volume of halite at the pressure of interest (Hu et al. 2021): The apparent partial molar volumes are a function of the composition of the solution, but also temperature and pressure (Appelo et al. 2014).We consider in the context of pressure solution, the molar volume of the halite V Halite to be constant at 27.1 cm 3 mol −1 , and leave the apparent partial molar volumes of Na + and Cl − to be calculated.The appar- ent partial molar volumes of the ions can be calculated as a function of the thermodynamic pressure and temperature using the formulation in Appelo et al. (2014), and this has been implemented into Crunch for the purposes of this study. (1) The kinetics of halite is a function of solubility K P,T , which is a function of stress and temperature.Using a transition state theory (TST) rate law with a rate constant, k , of 0.254 mol ⋅ m −2 s −1 from Alkattan et al. (1997), we obtain: where Q is the ion activity product defined by: The activity of the sodium ion, a Na + , is defined by the product of the activity coefficient and the concentration: and the activity of chloride, a Cl − , is defined similarly.
The rate law in Eq. ( 3) is combined with the diffusive flux described by Fick's Law, and an accumulation term to yield the reactive transport equations for Na + and Cl − linked together by the reaction rate of halite: where m D Na + and m D Cl − are the effective diffusivities (Steefel et al. 2015a) which are updated due to the change of the porosity as a result of deformation, grain movement, and chemical reaction (Hu et al. 2021).The symbol m is the Archie's Law cementation exponent that is introduced to capture the increased tortuosity at low porosity values (Steefel et al. 2015a;Hu et al., 2021;Steefel and Hu 2022).Here individual grid cells can have porosity values between 1.0 (representing the fluid phase) and 0.0 (pure solid) to represent the microscale structure (Li et al. 2008;Steefel et al. 2015b).Note that when the fluid phase with a porosity = 1 is considered, the cementation exponent has no effect, but when the porosity is low, then the cementation exponent results in a high tortuosity and low effective diffusivity.We use the approach of calculating diffusion in the solid halite grains rather than adopting complex no-flux boundary conditions along the grain boundaries.Assuming binary diffusion due to electroneutrality, a single value for both the diffusivity of Na and Cl is used.
Examining Eq. (3) in combination with Eqs.(1, 2), we see that the spatially variable solubility K P,T is higher when the compressive stress is higher.The resulting chemical potential gradients that develop drive diffusion of Na + and Cl − from highly stressed to relatively lower stressed contacts. ( The temperature and pressure dependence of the partial molal volumes of the sodium and chloride ions are calculated from the formulation in Appelo et al. (2014) reproduced in Hu et al. (2021).The temperature dependence of the halite solubility between 25 and 30 °C at standard pressure (1 bar) is calculated from the Livermore thermodynamic database (Zimmer et al. 2016).The temperature gradient along with the temperature dependence of the solubility results in a spatially variable chemical potential.Rates of reaction increase as well according to the activation energy for mineral dissolution.

Coupling Geometry
The challenges of linking of NMM to CrunchFlow are primarily associated with differences in the meshes and numerical interpolations in space and time.NMM uses a Lagrangian mesh and allows for boundaries across meshes, while CrunchFlow uses a Eulerian mesh that needs to conform to boundaries.Because of these differences, it is challenging to capture large deformation and moving boundaries due to contacts in both NMM and CrunchFlow, and transfer coupling terms between these two codes with sufficient resolution in space and time.
To overcome these challenges, we first developed a new approach in CrunchFlow to generate mesh automatically.With this approach, when the grain boundaries move due to deformation and/or motion, the Eulerian mesh in Crunch-Flow is reconstructed to capture such moving boundaries.In addition, to transfer values between NMM and CrunchFlow, we developed a mesh mapping algorithm to map between the rectangular CrunchFlow grids and the NMM irregular meshes, as shown in Fig. 3. Once we identify which NMM element (in yellow lines) contains a CrunchFlow cell center (blue nodes), we use a linear interpolation to calculate displacements at each cell center of CrunchFlow and associated stress tensor (i.e., a constant in each NMM element as a result of linear interpolation).To represent the arbitrary shapes of grain boundaries, we use a much denser mesh in CrunchFlow to be able to approximate each boundary of a grain.As shown in Fig. 3 (right), when the Crunch-Flow mesh is dense enough, we can approximate each grain boundary with a number of boundaries of different grid cells.Because we used a much finer mesh in Crunch and a coarser mesh in NMM, such approximation with a fine CrunchFlow mesh is sufficient to capture stress that is calculated by NMM even in CrunchFlow cells that are at boundaries.

Brine Migration in a Single Halite Grain due to a Temperature Gradient
As a first case we consider the migration of brine in a single halite grain in the presence of a temperature gradient.The conceptual model is shown in Fig. 4. In this example, For the given temperature gradient, we calculate the difference in temperature between the two sides of the inclusion.Using the thermodynamic model presented here and previously (Hu et al. 2021), we calculate the difference in solubility at either side of the inclusion.
The calculation domain is 4 × 2.4 mm with a 0.4 × 0.2 mm brine inclusion in the center.The left-hand side of the halite grain is 100 °C while the cooler side is 40 °C.The temperature difference across the fluid inclusion of 0.4 mm is thus 6 °C-this is an extreme case of temperature difference on the two sides of a fluid inclusion.The purpose of setting such a high temperature difference is to see how far the brine can migrate within a limited time.With a temperature of 100 °C at the high temperature side of the halite grain, an inclusion located in the center will have a temperature of 73.15 °C on the hotter side and 67.15 °C on the cooler side.For the solubility model considered here, this corresponds to a sodium concentration in equilibrium with halite of 6.5231 and 6.509 mol/kgw for the hotter and cooler sides of the inclusion respectively.The temperature dependence of the solubility between 40 and 100 °C is shown in Fig. 5.
The rate constant (Alkattan et al. 1997) used for halite dissolution and precipitation is quite high, so it appears likely that the rate-limiting process is diffusion across the fluid inclusion within the brine between (or close to) the equilibrium concentration values of 6.5231 and 6.509 mol/ kgw.(Fig. 6).The overall process of brine migration is slow since the entire halite lattice must be completely dissolved at the hot side and then reprecipitated on the cool end.Here we present simulations out to 10,000 h assuming an Archie's cementation exponent for calculating the diffusion tortuosity of 3.00 (Fig. 6) and 3.33 (Fig. 7).Theoretically, a higher cementation exponent of 3.33 leads to slower migration of the fluid inclusion up temperature.A distinct precipitation band at the trailing edge of the "bubble" migrating towards the hot end develops in the case where the cementation exponent is equal to 3.00 (Fig. 6), thus separating the original inclusion from the newly developed fluid inclusion.In both cases, the brine inclusion moves slowly up the temperature gradient since the entire lattice of the halite crystal must be dissolved.Note that the porosity range in both Figs.6 and 7 was selected to show the most intensive range of porosity changes, although the total range of porosity as a result of the dissolution is 1.0 to nearly 0.
Despite the fact that both the size of the brine inclusion relative to the size of domain and the temperature difference are extreme, the speed of brine inclusion migration is slow.Thus, we conclude that pure temperature may have a limited role in the dissolution and precipitation of a geometric feature such like a brine inclusion in a salt rock.Based on this, we will not show the detailed effects of temperature gradient on the migration of brine inclusion in the following section because these would be negligible for the temperature gradient and the sizes of the inclusions that are considered in those cases.

Effects of Temperature and Geometry on Pressure Solution in a Salt Rock
In a second set of simulations, we consider a grain pack consisting of halite derived from an image by Desbois et al.
(2010, Fig. 6a).In contrast with the case considered in Hu et al. (2021) where the grains including many sharp corners (i.e., first-order geometric features), the image presented here shows a salt rock with grains that have rounded and smoother asperities that are visible as second-order geometric features.The size of the domain is 7.5 × 2.6 mm.Based on a previous comparison between continuous and discontinuous approaches for simulating a fractured rock (Hu and Rutqvist 2022), we chose the discontinuous approach so as to be able to capture rigorously the contact stress distribution.
The geometric representation and initial meshes that are used in NMM and CrunchFlow are shown in Fig. 8.In both NMM and CrunchFlow simulations, we rigorously represent the second-order geometric features including the asperities along each grain boundary as well as the major holes or fluid inclusions within the grains.We assume all the minerals are 100% halite with a Young's modulus 3GPa and a Poisson's ratio 0.3.The friction angle of all solid surfaces is 30°.We apply a load of 5 MPa on the right boundary and a load of 3 MPa on the top boundary, respectively.The left and bottom boundaries are fixed.The cementation exponent for this set of simulations is assumed to be 4.0-this value ensures that there is relatively little diffusion into the solid halite grains.Because the system is tightly packed with small pore space between the grains, we do not observe large displacements in the system.After reaching mechanical equilibrium, the volumetric strain is shown in Fig. 9 and the distribution of minimum principal stress is shown in Fig. 10.We see that the second-order features dominate where contacts occur between every two surfaces of two grains.As a result, only the elements that are adjacent to the second-order features are highly compressed (shown in light to dark blue in Fig. 9).The stress concentration (shown in Fig. 10, light to dark blue) occurs not only in areas where contacts occur, but also at some other secondorder geometric features that are not in contact as well as at some of the major fluid inclusions.
To interpret the stress concentration in a grain at the bottom (we name this as "Grain A"), we show the zoomed-in geometry in Fig. 11.Note that all the grains at the bottom are truncated from a larger image by Desbois et al. (2010), (Fig. 6a).As a result of such truncation, Grain A has a relatively large area of geometric features within a small area of itself.These geometric features include two major brine inclusions, two geometric abnormal features, and relatively rough boundaries for such a thin and long grain (shown in Fig. 11).As a result of these geometric features, we see high stress concentration within the entire grain.It is important to note that such a grain geometry and stress distribution may be unrealistic due to the truncation.But here we make use of the magnified stress distribution and investigate how these intense geometric features and stress concentration affect dissolution, precipitation, and pressure solution.
For the purpose of understanding the relative importance of temperature and stress on pressure solution, we carry out simulations with the stress distribution shown in Fig. 10 and a linear temperature gradient from 25.0 to 30.5 °C (from right to left) over the 7.5 mm width of the domain.The pressure and temperature dependence of the halite solubility for these ranges are shown in Fig. 12.Note that for the pressure and temperature ranges considered, the change in halite solubility is very similar.The stress range is relatively small because the asperities (i.e., the second-order geometric features) are evenly distributed across the system.Thus, the solubility change within this stress range is nearly linear, which is different from the case of compaction of a salt aggregate considered in Hu et al. 2021.The temperature gradient is high compared to what one would expect in a natural setting, but this choice allows us to examine the behavior for an approximately equal effect of the two parameters (temperature and pressure) on halite solubility.Although including a temperature dependence for the diffusion coefficient is straightforward, we did not do so in these simulations because of the very minor effect for such a small (absolute magnitude) temperature difference-the effect is dwarfed by the temperature effect on the solubility.
We show results of pressure solution and precipitation in the following scenarios: (Case A) only temperature and no stress, (Case B) only stress and no temperature, and (Case C) with both stress and temperature.

Case A
As a first case, we consider the salt rock with a temperature difference of 5.5 °C, with the left-hand side at 30.5 °C and the right-hand side at 25.0 °C.No stress is considered.
Figure 13 shows the changes in porosity over time.Overall, we do not see visible changes in terms of pore space increase or decrease across the system.Instead, the changes are minor and where they occur, they are spread across the entire system.To better understand the detailed processes, we show the time evolution of the dissolution rate and the precipitation rate in Figs. 14 and 15, respectively.We use circles in green to outline the pore spaces that are separated by asperities of grains that are in contact.It is clearly shown that the temperature gradient across the entire system take effects on these individual pore spaces.These individual pore spaces then function like isolated "islands" with dissolution at the hotter ends (Fig. 14) and precipitation at the cooler ends (Fig. 15).When it reaches 3000 h, we see that the dissolution occurs deeper toward certain grain surfaces.However, the temperature gradient is not sufficient to dissolve any asperities that isolate these pore space.Despite these localized changes, the temperature gradient is not large enough to have a substantial effect over the entire system.As discussed further below, this is largely because the temperature gradient exists over the entire width of the domain, with no "hot spots" of dissolution developed (in contrast to the stress cases considered below).Other simulations (not shown here) demonstrate that higher temperature gradients can accentuate these effects, eventually leading to some brine migration towards the hotter end.But these temperature gradients are geologically unrealistic, so they are not considered further here.

Case B
In Case B, we consider the stress distribution shown in Fig. 10 under isothermal conditions (25 °C).As in Hu et al. (2021), we approximate the minimum principal stress as a thermodynamic pressure (Hobbs and Ord 2016), with the local brine in contact with the stressed solid having the same pressure.With this approach, we can carry out a rigorous thermodynamic calculation of the solubility of the two coexisting phases (halite and brine) at the same temperature and pressure.
Figure 16 shows the time evolution of the porosity due to the stress effects, with the zone of maximum dissolution corresponding to the maximum of the minimum principal stress shown in Fig. 10.As analyzed in Fig. 11, because of the dense distribution of geometric features in Grain A, stress concentration occurs in this entire grain.As a result of stress concentration, porosity increase first occurs at the geometric features where Grain A is in contact with the two upper grains (1000 h), and then this porosity increase continues propagating to a larger area (3000, 5000, 7000, and 10,000 h).
We see that the dissolution first occurs at the stress concentration areas of geometric features.This occurs across the entire domain.But now we focus on analyzing the area of Grain A where most of the dissolution evolution occurs.The geometric features in the area of Grain A include the two contact areas between Grain A and the upper two grains, and the two inclusions where compressive stress is concentrated and high.As a result of pressure solution, these geometric features are dissolved (100 h), thus gradually connecting the isolated pore space and creating an open pore channel for the dissolved solid mass to diffuse.Subsequently, lower stress starts to take effect so the two localized areas of pressure solution (i.e., contact areas and inclusions within Grain A) break through (1000 h) and propagate to a much larger area along different grain boundaries where stress is relatively higher.
The precipitation that occurs is more widely dispersed in the system (Fig. 17) but is primarily located along grain boundaries and boundaries of inclusions with lower compressive stress at the early stage (100 h).Thus, precipitation, pressure solution and diffusion reach mass balance in localized areas.With an increased area of pressure solution, we see that the precipitation starts to occur near areas that are newly dissolved (1000, 3000, 5000 h) until a maximum is reached.At the latest stage, precipitation starts to concentrate in the area of Grain A where a significant amount of solid mass has been dissolved (10,000 h).It is important to note that here the rates are instantaneous rather than cumulative rates.

Case C
In the third case we consider the evolution of the salt rock due to both stress and temperature.Given the modest temperature-only effects shown in Fig. 13, we can expect that the results will not be greatly different from what we observe in Case B, but these subtle differences may become significant as longer time scales are considered.
Figure 18 shows the time evolution of porosity.Overall, the porosity shows similar change patterns as in Case B. After a closer comparison between Figs. 16 and 18, we see the temperature gradient does not fundamentally change the porosity distribution, but shifts the maximum in the overall porosity change to the left side-the hotter side.This become particularly obvious after 7000 h.
Figures 19 and 20 show the results of halite dissolution rate and precipitation rate over time, respectively.Comparing Figs.21 and 19, we see that dissolution rate become noticeably different after 5000 h.Starting at 5000 h, the area near Grain A shows greater dissolution on the left.Starting from about 7000 h, the two contact areas on the left near the hotter side of the domain continue dissolving.Based on Figs. 18,19,20, we interpret that this is because pressure solution removes the solid mass of the contact geometric features, with the result that several isolated pore spaces become connected, giving a larger space for diffusion and precipitation.Under a higher temperature, this larger area that maintains localized mass balance results in a higher rate of pressure solution (shown in Fig. 19) and a higher rate of precipitation (shown in Fig. 20).It is interesting to note that at 10,000 h, the precipitation rate reaches zero in much of the system-which explains the increase of porosity on the left side of the system.In contrast, near the area of Grain A the precipitation rate increases as a result of intense pressure solution in that area.

Discussion
In the second set of simulations, we extracted an image from a natural salt rock and conducted THMC modeling by deactivating mechanical and thermal components separately to investigate the individual roles of temperature and mechanics, and their interplay with the natural geometry where a large number of second-order geometric features (e.g., asperities along the grain boundaries, inclusions, abnormal geometric features) are represented."Grain A" and other grains in the bottom are truncated from an original image.Because of that, and in combination with the boundary effect, we have unusually intense geometric features within this long and thin grain (i.e., Grain A).Because of this unique distribution of geometric features, we saw high stress concentration in the entire grain, and thus pressure solution effects that dominate in this part of the system and propagate to other parts of the system.If we had the other truncated half of the image below the bottom, we might not be able to obtain the patterns of stress and porosity as shown in 4.2.Because of these extreme conditions, the overall changes of porosity and permeability in such a system as a result of pressure solution may be unrealistic for the time scales considered.However, it is useful to consider such an extreme case with intense geometric features to show the significant roles of concentration of geometric features that lead to system dissolution, breaking through different inclusions, and forming a large pore space.
In both sets of simulations, we applied large temperature gradients-we applied a 6 °C difference for a 0.4 mmlong inclusion, and a 5.5 °C temperature difference for a 7.5 mm-long domain.Both are extremely large gradients that may not be realistic in most geological settings.The purpose of having such large numbers is to investigate how temperature gradients affect dissolution, precipitation, and pressure solution under the most extreme conditions.For the natural salt rock simulations, the choice of this temperature gradient also allows us to consider the effect of the two variables, temperature and stress, when their effect on halite solubility is approximately equal.Based on the simulation results, we saw temperature has a very limited impact in both systems.However, temperature may take effect if new surfaces are created-for example, via thermal fracturing.In that case, the geometry and stress distribution can be significantly changed, and this will affect the distribution of stress-induced pressure solution.
Because there are no existing experiments or numerical models available for validating the model, we have carried out a step-by-step and component-by-component self-validation in previous studies.Previously, the mechanical model and chemical model have been individually validated.Comparison with experimental studies in the future would need to consider the case where the halite is not 100% sodium and chloride, that is, pure halite (as we have assumed in this study).One expects that an impure halite would show a different solubility, and perhaps different rates of reaction.
It is interesting to note that the extreme case of stress concentration due to the geometric features and truncation associated with Grain A may be useful for analyzing potential mechanical-chemical processes in an excavation of a salt cavern (for example, for hydrogen storage).Despite the fact that the scale (microns) that is considered in this study is much smaller than the scale (meters) of an engineering setting for a tunnel, the relative scale of the geometric features to the domain size of interest are similar.As a result of excavation, stress concentration may occur in some of the excavation surfaces, leading to pressure solution.If other geometric features (such as fractures) are present, the pressure solution area may propagate, leading to a larger dissolved area in the salt cavern.In addition, the dissolved species may react with the stored gas (taking hydrogen as an example).Thus, it is important to consider the coupling of excavation-induced geometric and mechanical changes with the reactive transport, and to consider the effects of gas storage in salt caverns for engineering applications.

Conclusions
Pressure solution, a mechanism that involves tight coupling between the geometry and thermal-hydro-mechanical-chemical (THMC) processes, plays an important role in diagenesis.A number of experimental, theoretical, and numerical studies have been conducted on different aspects of pressure solution.However, a quantitative understanding and prediction of how geometry affects pressure solution in sedimentary rocks has never been attempted.From the perspective of numerical modeling, the challenges to achieving such a goal are associated with computational geometry and its tight coupling with THMC processes.In this study, we tackled these challenges and made the first attempt to conduct microscale THMC modeling to understand and quantify the impacts of geometry and temperature on pressure solution, using natural salt rock as an example.This modeling capability is achieved by expanding a novel MC code that we developed previously (Hu et al. 2021), this study by including temperature gradients.
To highlight the important roles of geometric features at different stages of diagenesis (i.e., from grain aggregates to sedimentary rocks), we first established the concept of first-order, second-order and Nth-order geometric features.In our previous study (Hu et al. 2021), we found that the first-order geometric features (i.e., the sharp corners of grains) play key roles in grain relocation and pressure solution in the process of compaction of a grain aggregate.Based on that finding in combination with the new definition of Nth-order geometric features, we developed a hypothesis: Depending on the stages of diagenesis (i.e., from grain aggregates to sedimentary rocks), first-order, second-order, and Nth-order geometric features sequentially play key roles in dominating contact dynamics, contact stress, and pressure solution in a sedimentary system.
To test this hypothesis and answer the scientific question "How do the geometry and temperature affect the pressure solution in a natural sedimentary rock such as a salt rock?", we conducted two sets of simulations with different scenarios.
In the first set of simulations, we used simple geometry-a single rectangular halite grain with a single brine inclusion in the center.It is assumed that because of the temperature difference, the brine inclusion would dissolve at the hotter face and precipitate at the colder face, resulting in migration toward the heat source.However, for the given temperature difference and the given domain size, and for the two different sets of parameters, we observed that the speed of brine inclusion migration is slow.Thus, we concluded that pure temperature may have a limited role in the dissolution and precipitation of a geometric feature such like a brine inclusion in a salt rock.
In the second set of simulations, we extracted geometry from a natural salt rock image and conducted simulations with different cases: (A) only temperature and no stress, (B) only stress and no temperature, and (C) with both stress and temperature.We observed different phenomena for these different cases.
In Case A, it is shown that the temperature gradient across the entire system has an effect on the local pore spaces that are separated by contacting asperities.These individual pore spaces thus function like isolated "islands" with dissolution at hotter ends and precipitation at cooler ends.Because the temperature gradient is not sufficient to dissolve any asperities that isolate these pore spaces, these localized dissolution and precipitation occur uniformly in the system where pore spaces are present.
In Case B, intense geometric features (e.g., major asperities, inclusions) in one area undergo stress concentration, thus dominating pressure solution in that area.Gradually, the isolated pore spaces are connected due to pressure solution, creating an open pore channel for the dissolved solid mass to diffuse.Subsequently, lower stress starts to take effects so the two localized areas of pressure solution break through and propagate to a much larger area along different grain boundaries where stress is relatively higher.
In Case C, on top of Case B, pressure solution is spread out at contacting highly stressed geometric features that are close to the hotter side.Thus, we see increased porosity on the left side of the system.
Based on these simulations, we found that: • Second-order geometric features dominate pressure solution in a natural sedimentary rock because they dominate the stress distribution.• Pressure solution increases the connectivity in a system, thus potentially increasing dissolution over a large area where stress is sufficiently high.• Temperature gradients may have a limited impact on pressure solution, but can change the system behavior by changing the magnitudes of pressure solution in localized areas.
Thus, our microscale THMC model has tested the hypothesis regarding the roles of second-order geometric features on a sedimented salt rock.And this new model has provided quantitative descriptions of how geometry and temperature impact pressure solution.In the future, we will continue testing this hypothesis with more sets of samples in combination with field and laboratory observations.Based on the findings from this study, we will make use of geometric features as signatures for identifying the stages and key processes of the Earth's sedimentary rocks.We will also make use of the toolset that was developed for pressure solution to analyze MC processes in salt caverns for gas storage (such as hydrogen storage in salt caverns).

Fig. 7 Fig. 9
Fig. 7 Brine-filled fluid migration up a temperature gradient.Top: Halite rate in mol m −3 s −1 , with positive values representing precipitation, negative values dissolution.Bottom: Time evolution of the

Fig. 10
Fig. 10 Distribution of minimum principal stress in the system (Unit: Pa)

Fig. 12
Fig. 12 Halite solubility mol/kg water (molality).Left: Pressure dependence from 0 to 120 bars pressure-a relatively small range of stress resulting in a linear solubility distribution.Right: Temperature dependence between 25 to 30 °C

Fig. 16
Fig. 16 Time evolution of the porosity for Case B. Times from left to right and top to bottom are 100, 1000, 3000, 5000, 7000, and 10,000 h

Fig. 21
Fig. 21 Time evolution of the halite precipitation rate for Case C. Times from left to right and top to bottom are 100, 1000, 3000, 5000, 7000, and 10,000 h