Kinematics of Crystal Growth in Single-Seal Syntaxial Veins in Limestone - A Phase-Field Study

Building on recent developments in phase-field modeling of structural diagenesis, we present an analysis of single-seal syntaxial calcite vein microstructure in a variety of limestones. We focus on the effects of fracture aperture, intergranular versus transgranular fracturing, crystal habit and the presence of second phases in the host rock, to systematically investigate a simplified set of models covering the main classes of limestone in 2D. We incorporate the kinematic process of growth competition between differently oriented crystals, growth rate anisotropy between rough and faceted crystal surfaces and different growth rates on intergranular to transgranular fractures. Results show that within the considered parameter space we can reproduce a wide range of vein microstructures in limestone known in nature, such as stretched crystals, wide-blocky veins, and elongated crystals. We identify five archetypes of vein microstructures in limestones, which are diagnostic for different kinematics and evolution of transport processes and illustrate the effect of key parameters in microstructure maps. We show how syntaxial veins with median line form after intergranular fracturing, while stretched crystals indicate transgranular fracturing. Intergranular fracturing leads to stronger growth competition and more prominent CPO in syntaxial veins. Our results can


Introduction
The chemistry and mechanics of fracture sealing in deforming hydrothermal and reactive environments is important in basic geoscience and applied studies. Examples are the seismic cycle (e.g., Alexandrakis et al., 2014), evolution of ore deposits (e.g., Cox et al., 1987;Sibson et al., 1988), geothermal energy (e.g., McNamara et al., 2016), and the production and injection of fluids from and into the subsurface (e.g., Almansour et al., 2020).
In dilatant fractures, aqueous fluids can dissolve or precipitate minerals, and this in turn affects the transport and mechanical properties of the rock system. Syntectonic syntaxial calcite or quartz veins are ubiquitous in rocks at low metamorphic grades in the upper and middle crust (Passchier & Trouw, 2005;Ramsay & Huber, 1983). The formation of these veins includes multiple coupled mechanical, hydraulic, Abstract Building on recent developments in phase-field modeling of structural diagenesis, we present an analysis of single-seal syntaxial calcite vein microstructure in a variety of limestones. We focus on the effects of fracture aperture, intergranular versus transgranular fracturing, crystal habit and the presence of second phases in the host rock, to systematically investigate a simplified set of models covering the main classes of limestone in 2D. We incorporate the kinematic process of growth competition between differently oriented crystals, growth rate anisotropy between rough and faceted crystal surfaces and different growth rates on intergranular to transgranular fractures. Results show that within the considered parameter space we can reproduce a wide range of vein microstructures in limestone known in nature, such as stretched crystals, wide-blocky veins, and elongated crystals. We identify five archetypes of vein microstructures in limestones, which are diagnostic for different kinematics and evolution of transport processes and illustrate the effect of key parameters in microstructure maps. We show how syntaxial veins with median line form after intergranular fracturing, while stretched crystals indicate transgranular fracturing. Intergranular fracturing leads to stronger growth competition and more prominent CPO in syntaxial veins. Our results can be extended to 3D to include multiple crack-seal events, pore-space cementation and simulation of fluid flow, providing a generic platform for modeling structural diagenesis in limestones.
2 of 23 and geochemical processes and parameters (e.g., temperature, chemical composition, pressure of fluid) and results in a broad spectrum of macro-and microstructures (Laubach et al., 2019). Veins can be used to reconstruct the history of stress and fluid pressure (Bons et al., 2012). Chemical and isotope analysis of vein cement (Arndt et al., 2014) and fluid inclusions provide evidence on the thermal and chemical environment during vein formation and can be used to draw conclusions on fluid flow and mass transfer in the system (Boullier & Robert, 1992;Cox, 2007;Fisher & Brantley, 1992). However, the evolution of vein microstructures (Durney, 1973;Passchier & Trouw, 2005;Ramsay & Huber, 1983) remains an open field (e.g., influence of fracture type). Urai et al. (1991) presented the first model to explain the kinematics of crystal growth into irregular spaces and the evolution from euhedral into nonSchmidegg shapes (Schmidegg, 1928). This work was extended by Nollet et al. (2005), who explored the transition from crystallographically rough surfaces to faceted crystals and back to irregular shapes, when crystals grow against a rough wall-rock surface. Depending on growth direction (syntaxial, antitaxial, bitaxial, unitaxial Durney, 1973), fracture location (localized, delocalized), opening trajectory of the fracture, fracture type (intergranular, transgranular), and mineral type (e.g., calcite, quartz) a broad variety of vein morphologies with different characteristic microstructures can form (e.g., fibrous, elongated-blocky, stretched, ataxial crystals) (Bons et al., 2012). Becker et al. (2011) demonstrated how the presence of a second mineral affects the crystal morphology in syntaxial bi-mineralic multi crack-seal veins.
In syntaxial veins crystals grow from both sides of the fracture wall inwards into the center, whereas in antitaxial veins crystals grow epitaxially on one side of the fracture from a median zone in the vein center outwards toward the inert wall rock surface (Durney, 1973) (Figure 1). Since antitaxial veins have a different mineralogy than the wall rock (lower adhesion), the fracture opens at the vein-host rock interface and the opening and sealing occur continuously (e.g., due to force of crystallization) or in a sequential pattern. When opening increments are small, the development of crystal facets is limited and fibrous crystals can form (Ankit et al., 2013;Ankit, Selzer, et al., 2015;Bons et al., 2012;Hilgers et al., 2001;Nollet et al., 2005;Urai et al., 1991). The numerical modeling of antitaxial veins has two big advantages: (1) the sealed crack is very weak and will reopen in the same location so that the fracture mechanics is simple, and (2) crystals grow against an inert wall, so that the joining of crystals growing from both sides does not have to be simulated. Urai et al. (1991) developed a kinematic model for the formation of these veins and investigated the effect of the fracture wall roughness on the evolving microstructure. Expanding on this work, sharp interface methods, which explicitly track the growth front, have been utilized in 2D to understand the influence of opening velocity and trajectory in antitaxial veins (e.g., Vein Growth (Bons, 2001;Hilgers et al., 2001) or FACET (Zhang & Adams, 2002)).
The vast majority of veins in limestone consist of calcite which grows epitaxially on the calcite grains in the host rock. They are called syntaxial or ataxial veins due to the crystal growth from both sides toward the fracture interior (Hilgers et al., 2001;Passchier & Trouw, 2005). For simulations, this has two additional problems: (1) the sealed crack is strong and can reopen at a different location so that mechanics of fracturing is complicated and can be delocalized, and (2) crystals grow from opposite sides of the crack, so that simulations have to be performed against dynamic interfaces, instead of inert fracture walls.
The cellular-automaton-type program Prism2D (Lander et al., 2008) has been applied for crystal growth in syntaxial veins (in quartz: Lander et al., 2008;Lander & Laubach, 2015, in dolostone: Gale et al., 2010. In their numerical model they incorporated the well-known effect that crystals, which grow on rough surfaces, show a significant decrease in the total growth velocity after they attain their faceted shapes (see e.g., Periodic Bond Chain theory from Hartman & Perdok, 1955a, 1955b, 1955c. Their numerical model has been calibrated with quartz growth kinetics from laboratory experiments (Lander et al., 2008).
The phase-field method is widely used in material science for modeling microstructure evolution processes. This method uses a mathematical model to represent an interface, in this case for growing crystals (Fix, 1982;Langer, 1980Langer, , 1986. It is based on the local energy functional minimization and allows with its diffuse interface approach an elegant treatment of moving boundary problems where no explicit interface-tracking algorithm is required. Recently, the phase-field method has emerged as a powerful tool to model cementation of fracture and matrix porosity in the subsurface (see e.g., Ankit, Urai, & Nestler, 2015;Hubert et al., 2009;Wendler et al., 2016). It allows the investigation of different vein mineralogies, such as potash alum (Ankit et al., 2013), quartz (Ankit, Urai, & Nestler, 2015;Prajapati, Abad Gonzalez, et al., 2020;Wendler et al., 2016) or calcite (Prajapati et al., 2018;Spruženiece et al., 2020Spruženiece et al., , 2021 in 2D and 3D. In one of the first phase-field studies of antitaxial vein formation Hubert et al. (2009) utilized a nonfaceted smooth anisotropy formulation in the interfacial energy term for the modeling of antitaxial vein growth, whereas in nature crystals are flat faceted crystals with sharp corners. This limitation was overcome by the work of Ankit et al. (2013), who utilized a facet-type anisotropy formulation in the interfacial energy term (see e.g., Nestler et al., 2005). In this and a subsequent work (Ankit et al., 2013;Ankit, Selzer, et al., 2015) studies of antitaxial vein formation were conducted in 2D and 3D with varying fracture roughness, opening velocity, increment and trajectory, and showed a sound agreement with previous works of Urai et al. (1991) and Hilgers et al. (2001). However, antitaxial veins are relatively rare in limestones. They typically form in shales, marls and silty mudstones and consist of calcite, gypsum and halite (e.g., Cobbold et al., 2013;Hilgers & Urai, 2002;Leitner et al., 2014;Meng et al., 2019;Philipp, 2008).
The same phase-field model as in Ankit et al. (2013) has been applied in Ankit, Urai, and Nestler (2015) for studying in 3D multi-crack-seal events in syntaxial quartz veins with normal and oblique openings, and in Prajapati et al. (2018) for modeling calcite vein formation. Therein, the evolution of different crystal morphologies with varying fracture opening velocities is demonstrated. In Wendler et al. (2016) the phase-field model parameters were calibrated to laboratory quartz growth experiments (e.g., Okamoto et al., 2010) by using anisotropy formulations in the interfacial energy term and the kinetic mobility. In order to incorporate the effect of different growth velocities for rough and faceted crystals the kinetic mobility parameter in the phase-field evolution equation was adapted. Additionally, they utilize an interpolation of the mobility in solid-liquid multiphase regions to avoid artifacts by triple point pinning. However, this approach is limited to grain structures, where all crystals have similar total growth velocities.
Recent studies of syntaxial calcite veins in limestones by Spruženiece et al. (2020Spruženiece et al. ( , 2021 showed that in some veins the crystal growth rates vary depending on how the host rock grains are fractured (intergranular or transgranular). They utilized the anisotropy formulation of the kinetic mobility of Wendler et al. (2016) and combined it with the evolution equation of Steinbach (2009) to distinguish between fast-growing calcite crystals on transgranular fracture surfaces and slow-growing calcite crystals on intergranular fractures.
Limestones can have different mineral compositions and a wide spectrum in grain size distribution and crystal morphology. The grain-and crystal size distribution can vary from mm range (e.g., in crystalline limestone) to a few E m in e.g., mudstone (Dunham, 1962). Depending on the distribution of skeletal fragments (e.g., foraminifera, coral) the grain size differs from approximately equally sized grain structures (e.g., micrite, ooids) to mixed-sized structures, in which both fine and coarse grains occur (e.g., a fine-grained mud-dominated, nongrain supported limestone with at least 10% skeletal components, thus a wackestone). Figure 2 illustrates a diversity of veins in different limestone host rocks: mudstone, wackestone, and crystalline limestone. The vein microstructures vary from homogeneously sized crystals (stretched or bitaxial) over wide blocky crystals to elongated bladed crystals, depending on the host rock type and vein apertures.

of 23
In addition, depending on the host rock composition and fracturing conditions, the fracture surface can be highly heterogeneous, exposing varying amounts of intergranular and transgranular segments (Figure 3), as well as surfaces of accessory minerals, which are unfavorable for the nucleation and growth of vein crystals (Aubert et al., 2020;Spruženiece et al., 2020).
Even though previous studies examined how the fracture surface characteristics affect crystal growth in veins (Ajdukiewicz & Larese, 2012;Lander & Laubach, 2015;Lander et al., 2008;Urai et al., 1991;Williams et al., 2015), they mostly focused either on the fracture geometry or growth kinetics of single crystals.  Figure B1 in Arndt, 2016). wbc: wide-blocky crystal; hsbc: homogeneously sized blocky crystal; ebc: elongated bladed crystal; sc: stretched crystal; gc: growth competition. In this work, we build on these previous works aiming to provide a systematic framework for understanding the diversity of vein microstructures in limestone host rocks ( Figure 2) and generalize the occurring morphologies. We incorporate diverse limestone host rock microstructures and complex fracture cements in our models, similar to those that occur in natural limestones. We focus on the formation of syntaxial calcite veins in limestone and investigate single crack-seal events, where, after the initial opening and sealing, no repeated fracturing occurs.

Multiphase-Field Model
A thermodynamically consistent multiphase-field model (Nestler et al., 2005) is utilized to computationally investigate, characterize and classify crystal growth in limestone veins. The presented model was applied in the previous works of Spruženiece et al. (2020Spruženiece et al. ( , 2021 and is briefly described in this section. For further discussions on the model equations interested readers are referred to for example Ankit, Urai, and Nestler (2015); Wendler et al. (2016); Prajapati, Abad Gonzalez, et al. (2020).
The phase-field model is based on local minimization of the Helmholtz free energy, which is given by ( ) The temporal evolution of each phase-field E   is described with the variational derivative of the energy functional (Equation 1). We adopt the approach of Steinbach (2009) to apply different mobilities (growth velocities) for different phases without interpolation difficulties of the kinetic mobility in multi-phase regions. The evolution equation reads as with the kinetic coefficient 0 E M  and the kinetic anisotropy kin ( ) E a  n . For modeling a faster growth rate of rough crystal surfaces and a drop of the growth velocity after crystal facets have formed (Lander et al., 2008) we follow the approach of Wendler et al. (2016) and utilize for the kinetic anisotropy. The terms max  denote the largest and the second largest argument of the scalar products. When the crystal is in its equilibrium shape, max  are equal and the second argument (after E  ) vanishes, whereas when the crystal is in a nonequilibrium/rough shape the largest and second largest argument differ. The drop of the growth rate between rough to fully faceted can be calibrated to experimental data with the factor E  (see e.g., Prajapati, Abad Gonzalez, et al., 2020; Wendler et al., 2016 for quartz).

Computational Treatment With Assumptions
The model equations are implemented in the in-house software package Pace3D (Parallel Algorithms for Crystal Evolution in 3D), which is written in language C. Interested readers are referred to Hötzer et al. (2018) for a detailed description of implemented optimization and parallelization algorithms.
The model parameters used in the simulations are given in Table 1. We choose the simulation parameters in order to ensure a numerically stable simulation, while reducing computational costs. When additional physical properties (e.g., chemical composition, pressure, temperature of the fluid) are available the nondimensional simulation parameters could be mapped to physical values and information about the growth kinetics could be obtained (see e.g., Prajapati, Abad Gonzalez, et al., 2020;Wendler et al., 2016).
We assume, that precipitation is the rate-limiting step for the crystal growth, i.e., the attachment kinetics of the solute are slow compared to transport rate. This results in a constant supersaturation in the whole domain. Therefore, we utilize a constant chemical driving force bulk E f  in the present work.

Modeling Crystal Shapes
The major minerals in limestone are calcite, dolomite, and aragonite, and these can have a wide range of crystal habits (Goldschmidt, 1913). In the present work we focus on calcite that is the most common component in limestone. Calcite can develop different crystal habits, depending on temperature, supersaturation and impurities in the solution, which can influence the vein microstructure (Gutjahr et al., 1996;Nollet et al., 2005;Paquette & Reeder, 1990;Plummer et al., 1979). For this study we choose rhombohedral and scalenohedral crystal habits, as they commonly appear in nature (Bishop et al., 1999;White & Culver, 2011).
For the rhombohedral habit we apply the vectors in Equations 4 and 7 from Prajapati et al. (2017) with a c/a ratio of 0.86 (also given in Supporting Information S1), whereas we elongate the c-axis for the scalenohedral shape to obtain a c/a ratio of 3.41 in the 3D shape (e.g., Barber & Wenk, 1979;Gunasekaran et al., 2006;Sekkal & Zaoui, 2013). In Figure 4b the simulated 3D shapes for scalenohedral and rhombohedral calcite are depicted with the projections used in the 2D simulations (in a)). In the 2D projections the c/a' ratio is 1.02 and 3.32 for the rhombohedral and scalenohedral shape respectively (in Figure 4c).

Two Dimensional vs. Three Dimensional Models
Most of the simulations in the present work are performed in 2D to limit computational costs for setups with a large number of grains (from 890 to 14,660 grains). The 2D simulations are able to capture most aspects of the microstructure evolution, however the crystallographic orientation of the grains is only rotated perpendicular to the simulation plane and the third growth direction (out of plane) is not considered.
Here we note that at higher computational cost, an extension to 3D in general is possible as presented in Section 3.6 in a representative setup.

Parameter Space and Boundary Conditions
In this study we investigate the influence of the following parameters on the grain structure in the veins: aperture, fracture type, grain size distribution, crystal habit, and mineralogical composition of the host rock. Table 2 illustrates the combinations of the parameters explored in the following subsections.
In order to easily compare the simulated veins to natural limestones, in which the host rock grain size can have different scales (crystalline limestone in mm range; mudstone in E m range (Dunham, 1962)  Aperture The initial grain structures in the 2D and 3D computational domains were generated in a three-step process: (1) A host rock was generated with a Voronoi algorithm, in which the grain orientations are randomly distributed. Optionally, differently sized grain structures were combined to obtain the final host rock. (2) A fracture was inserted in the host rock, which cuts the grains either intergranularly or transgranularly. The fracture was created with a graphics editor and then converted into the data format of Pace3D. (3) The lower fracture wall was then shifted downwards in a direction perpendicular to the wall until the designated value for the opening aperture (e.g., 8 m D E ). The opening of the fracture was performed in a single step and the crystal growth was initiated after the opening of the fracture, which corresponds to a faster opening than crystal growth rate. The phase-field models are set up such that the simulations stop, after the vein is completely sealed.

Effect of the Fracture Surface Microstructure
In this section, we analyze different types of fracture surfaces that can occur in limestones, formed during intergranular, transgranular, and mixed type microcracking ( Figure 3). We assume a growth rate difference between intergranular and transgranular fractured grains, which can result from clay coatings on intergranularly fractured grains, in which a fracture cuts through the clay and leaves a full or partial coating on the calcite grains. This remant (partial) coating hinders the precipitation rate compared to freshly broken transgranular fractures. In contrast, a growth rate difference between different fracture types is not present, when, for example, a fracture cuts in between a grain and its coating or when the rock went through considerable primary pore space cementation and previous coatings are overgrown. We model the reduced growth rate on intergranular fracture surfaces as a decrease in the solid/liquid interface mobility 0 describes the ratio between the solid/liquid mobility of transgranularly to intergranularly fractured grains (abbreviated with TFGs and IFGs in the next sections). A higher E  corresponds to a larger growth rate difference, whereas E  = 1 implies no growth rate difference. In the forthcoming sections we set E  = 5 if not stated otherwise, which gives a sound agreement of the models to natural veins in micrite (Spruženiece et al., 2020).
In the next sections we refer to crystallographically controlled growth competition, when favorably oriented crystals outgrow and terminate unfavorably oriented crystals, whereas fracture-surface controlled growth competition occurs, when E  1 and slow growing IFGs are terminated by faster growing TFGs.
The simulation setup (160 × 2500 E x) for the investigation of the fracture type is depicted in Figure 5a. The percentage of intergranular fractured grains is varied from 0% to 100% and apertures of 2, 8, and 30 m D E are examined. Furthermore, a computational fracture with an aperture of 300 m D E with (1) 100% intergranular fracturing, and (2) mixed type fracturing ( E 87% IFGs and 13% TFGs) is generated (10000 × 4600 E x with E  6,300 grains in Figure 7; Setup depicted in Supporting Information S1).
In the next sections we distinguish between the following crystal shapes: stretched-equidimensional (crystal height (h) approx. crystal width (w)), stretched-elongated (h E  w), (stretched) wide-blocky (w E  h), (bitaxial) blocky (h approx. w), and (bitaxial) elongated blocky (w E  h). Stretched indicates that two fragments of a broken crystal from the opposite fracture wall join, whereas bitaxial indicates that different crystals reconnect in the vein center. Note: The term stretched is also used in multi-crack seal veins for crystals which are fractured and reconnect multiple times, whereas in this work only one sealing event occurs.

Pure Intergranular vs. Pure Transgranular Fractures
First we compare the crystal morphologies in a computational fracture with apertures of 2, 8, and 30 m D E , where either pure intergranular or pure transgranular fracturing occurs (Figures 5b and 5c). When all grains are fractured along the grain boundaries (100% intergranularly), the vein crystals grow euhedrally on the surfaces of the wall rock calcite grains and a distinct median line forms for all tested fracture apertures (highlighted with black line) and stretched crystals are absent. For small apertures (2 m D E ) a majority of crystals reach the vein center, where they terminate against differently oriented crystals from the opposite side of the fracture wall. When the aperture increases, more crystals are terminated at early stages of the vein 9 of 23 evolution as they are outcompeted by their favorably oriented neighbor crystals. In veins with the largest aperture, crystals with elongated and blade-like morphologies dominate but there is still a median line.
In the veins with pure TFGs, most crystals tend to bridge the vein in fractures with the smallest aperture (2 m D E ), as the transgranularly broken wall rock grains have the same crystallographic orientation on both sides of the fracture, resulting in stretched crystal morphologies. With an increased initial fracture aperture (e.g., 8 m D E ) crystallographically controlled growth competition dominates and unfavorably oriented crystals are terminated. In some areas, when nonidentical crystal fragments meet in the vein center, a poorly defined median line is suggested. With a larger aperture more nonidentical fragments touch and the median line becomes more prominent. In the aperture of 30 m D E still a few stretched crystals form and this structure shows similarities to the pure intergranular fractured veins.

Mixed Fracturing Type
In this section we analyze the simulations for mixed type fractures with varying ratios of intergranular and transgranular fracture segments. Based on previous studies  we assume different growth rates for the two types of fracture surfaces, which results in fracture surface-controlled growth competition. Representative cut-outs of the fully sealed veins are depicted in a morphology diagram in Figure 6a and the mean crystal size is plotted in Figure 6b. Therefore, the vertical simulation domain width (scaled by m D E ) is divided by the number of crystals at the median line. The dominant crystal shapes in veins are classified in Figure 6c (see definitions in previous section). The complete images of the computational veins are given in the Supporting Information S1. Figure 6d shows the numbers of surviving crystals against the distance from the vein center (raw data given in repository). (Here, the fully sealed vein is equidistantly sliced and the number of occurring crystals in each slice is counted). , where crystallographically controlled growth competition between different TFGs and early termination of IFGs occurs. This leads to only few stretched-elongated crystals and more elongated blocky crystals.
In the veins with the aperture of 300 m D E (Figure 7) similar mechanisms occur. In the simulations on mixed fracturing type, IFGs are terminated at early stages and further crystallographically controlled growth competition proceeds exclusively between TFGs. The greater the distance between IFGs and TFGs is, the more The grain orientation in the host rock and along profile lines 1-4 in the vein is depicted as a rose diagram in Figure 7c. In mixed-type fractures, due to their increased growth rate, TFGs expand laterally and survive even if they are unfavorably oriented, which also results in an inhomogeneous growth front. In contrast, when all grains fracture intergranularly, only favorably oriented crystals survive and terminate less favorably oriented crystals close to the rim (Figures 7a and 7c). Here, the lattice preferred orientation in the vein is much stronger (profile line 2). The surviving crystals form elongated-blocky and fibrous crystals and are narrower than in the vein with the mixed type fracturing. The growth front in veins on pure IFGs tends to stay more regular and homogeneous.

Effect of Growth Rate Ratio   Between TFGs and IFGs
In this section, we vary the relative growth rate of TFGs with respect to IFGs by assigning a coefficient E  = 1, 5, 10, 20, inf., where E  = 1 implies no coating of IFGs. The term inf. indicates completely coated intergranular fractured calcite grains, which do not grow. The steps in between describe partly coated intergranular fractured calcite grains. We use a simulation setup with fracture apertures of 2, 8, 18, and 30 m D E . (220 × 2500 E x), where the black dots indicate the location of the TFGs. Additionally, we highlight accessory minerals (e.g., quartz grains) in the simulation setup in light gray in the forthcoming sections. Since the fluid in the open fracture is assumed to be supersaturated with respect to calcium carbonate and the vein crystals can only grow on calcite host rock grains, these accessory minerals are treated as inert (non growing) phases. m, as in mudstone, up to several mm as in crystalline limestone (Dunham, 1962;Flügel, 2013). In addition, many of the limestone types consist of a mix of fine and coarse grains, as in wackestone or packstone, where various proportions of mm-sized grains are embedded in a finer-grained matrix.
To test the effects of the host rock grain size on the vein microstructure, we compare two types of host rock. One case is set up with similarly sized grains, representative of crystalline limestone or mudstone structure in limestones and the other case, with differently sized grains representing a wackestone structure (Figure 9). The boundary conditions of these simulations follow those as in Figure 5a. To reduce the computational costs a smaller domain (12450 × 900 E x) was cut out (Figure 9a), in which approx. 8% of the volume comprises coarse grains (total of 14,460 grains in domain). Furthermore, a setup with a single big grain was extracted from the structure for the simulation with small apertures to further reduce computational time.
For both host rock types we investigate two cases, 100% transgranular fracturing and mixed type fracturing where only coarse grains are fractured transgranularly. In the veins with 100% transgranular fracturing (Figure 9b) similar vein morphologies develop for both types of host rock as in Figure 5c, namely stretched crystals in small apertures and bitaxial crystals in larger apertures. As expected, because of the absence of other grain boundaries, the coarse grain (in center) survives in all of the tested apertures, while the crystal width does not change significantly.
In the veins with mixed type fracturing ( E  = 5) the crystal morphologies that develop show the same trend of an increased termination of IFGs with increasing aperture in both host rock types as in the previous sections ( Figure 6). However, with larger apertures ( E 8 m D E ) the crystal morphologies in the veins differ in the two host rock types with the same aperture (Figure 9c), since the distance between TFGs is larger for the different sized grained host rock.

Effect of Inert Second Phases
Most limestones are not monomineralic. They can contain up to 50% of accessory phases, such as clay, quartz grains, feldspar and organic material. As the setup for all simulations assumes that the fluid in the fractures is supersaturated with respect to calcium carbonate, the vein crystals grow exclusively on the host SPÄTH ET AL.

10.1029/2021JB022106
14 of 23 rock calcite grains that are exposed on the fracture surface. Therefore, we set all accessory minerals to be inert for crystal growth and mark them in the simulation results in light gray.
A morphology map with representative cut-outs from the fully sealed veins is depicted in Figures 10a  and 10b for 82% and 100% intergranular fracturing respectively. For varying percentage of accessory minerals, we observe similar crystal morphologies as in the previous sections. However, in the small veins (2 m D E ) the crystals become wider, whereas for larger apertures the crystal morphology in the vein interior is not significantly impacted by the amount of accessory minerals (Figures 10a-10c). Overall, the effect of inert/ accessory minerals seems to be small for the tested cases, where they constitute up to 50% of the host rock grains.

Effect of Crystal Habit
Calcite in nature has a variety of crystal habits, which is controlled in a complex way by growth conditions, supersaturation and fluid chemistry. In order to determine how the microstructural evolution of calcite veins depends on the specific crystal habit, we investigate the crystal growth in rhombohedral and scalenohedral habits, which commonly occur in nature (Bishop et al., 1999;White & Culver, 2011). In addition, a variety of fracture apertures and fracture surface types are simulated with these crystal habits. Figure 11a depicts representative cut-outs of rhombohedral and scalenohedral calcite morphologies on mixed type fractures (82% intergranular) using the same initial setup as for the simulations in Figure 8a. Overall, similar wide-blocky microstructures form for both crystal habits. The main differences are caused by the more elongated (fast-growing) c-axis for the scalenohedral habit. As a result, the unfavorably oriented crystals are less likely to bridge the vein at large apertures, even if they grow on TFGs (Figure 11a ). Thus, in the vein microstructure for the scalenohedral case, the size and shape of the bridging crystals is less homogeneous, compared to the rhombohedral cases. Here, slow-growing grains on IFGs bridge the veins only for the small fracture apertures (Figure 11a) ). However, in simu-   Figure 11c, where an increased crystal size of the scalenohedral IFGs is observable. Figures 7a and 12a show simulations of rhombohedral and scalenohedral calcite morphologies on 100% intergranular fracture surfaces for large apertures. In both cases the final vein microstructure consists of highly elongated crystals (with aspect ratio E 10), that show crystallographically preferred orientation (Figures 7a and 7c-profile 2; Figure 12b) and a clear median line. In both cases the crystal size significantly increases toward the vein center, while the number of crystals decreases (Figure 12c), indicating crystallographically controlled growth competition, that is less pronounced on mixed type fracture surfaces (Figure 7b). Interestingly, for the scalenohedral habit the preferred orientation of the surviving large crystals is different compared to the simulations with rhombohedral habit (Figures 12b and 7c). This is attributed to the c/a ratio differences between the two habits used in the simulations. The final crystal structure looks similar for both habits, however the crystal shapes are blockier for the scalenohedral habit, whereas the crystals in the rhombohedral vein look more fibrous (Figures 7a and 12a). The growth front is less homogeneous, the median line is more serrated and the angles at triple grain-grain-liquid points during the grain growth are more acute in the scalenohedral case, compared to the simulations with the rhombohedral crystal habit.

Extension to 3D
An extension of the phase-field simulations from the 2D to 3D is straight forward, albeit increasing computational time. We perform an example simulation of a 3D vein in a host rock, which corresponds to a wackestone type in limestones. As in the 2D simulations (Figure 9), we assume, that the coarser host rock grains fracture transgranularly, whereas the finer-grained matrix fractures intergranularly. The fracture aperture is set to the size 32 m D E similar to the 2D simulation in Figure 9c and the growth rate ratio E  is set to 5. The generated domain (see Figure 13a) has the size of 2400 × 700 × 600 E x and comprises over 30,000 grains. Similarly to the 2D experiments, crystals that grow on TFGs are the first to bridge the vein, joining up with their fragment that grows from the opposite side of the fracture; and then expanding laterally, while terminating the neighboring IFGs that grow at a slower rate (Figure 13b). The fully sealed vein in 3D is depicted in Figure 13c and along the cut-outs and in Figure 13d. All crystals that bridge the vein have originated on the TFGs (illustrated with low opacity in c), whereas crystals, growing on IFGs form a rim of small crystals along the vein margins, where grain size increases toward the boundary between two TFG grains, similar to the observations in the 2D simulations (Figure 9c). In contrast to the 2D models, in 3D the crystals can grow in all directions. This can be seen in the 2D cutouts of the 3D simulation ( Figure 13d) as crystal "fragments" that are not connected to the host rock. Overall, however, this effect is small and does not significantly affect the microstructure of the vein.

Discussion
In this study, we systematically evaluated vein formation in a diverse range of limestones. Our numerical results show many similarities with natural vein microstructures and builds on the previous numerical works of Nollet et al. (2005); Lander and Laubach (2015); Prajapati et al. (2018); Spruženiece et al. (2020Spruženiece et al. ( , 2021, where controlling factors in vein formation are discussed. Based on the results of our simulations, we identify five distinct crystal archetypes: stretched equidimensional, stretched elongated, wide-blocky, (bitaxial) blocky, (bitaxial) elongated blocky crystals. The results suggest that the most important factors which determine the vein microstructure are the type of fracturing (IFGs and TFGs, when E  1) and the initial aperture of the fracture. In syntaxial veins stretched crystals form when transgranular fracturing occurs and the two parts of the same broken grain are sealed again by vein material, whereas bitaxial crystals grow either on intergranular fractures or misoriented transgranular fractures, where during the sealing process two unrelated crystals with different orientations join. The size of the initial fracture aperture determines if the vein is dominated by blocky or elongated crystals. Wide-blocky crystals form as stretched crystals that expand laterally due to the growth rate differences between intergranularly and transgranularly fractured grains (when E   1). However, the formation of wide-blocky crystals is limited to special cases of mixed type fractures, where only 5%-20% grains have a faster growth rate than the rest (medium aperture Figure 6a). With an increased aperture, elongated blocky crystals form instead of wide-blocky crystals, since after the early 18 of 23 termination of IFGs, normal growth competition takes place between the fast-growing crystals on TFG surfaces (Figure 7b).
We focus on syntaxial vein formation in fractures, where the opening velocity of the fracture is higher than the velocity of crystal growth, so that crystals can obtain facets. This is a common scenario for vein formation in nature. Blocky crystal morphologies and growth competition structures, such as gradually increasing crystals sizes toward the vein interior due to the termination of unfavorably oriented crystals are considered to be a strong evidence for veins that form in open fractures (Bons et al., 2012). Our results, however, indicate that crystal growth in open fractures can also result in stretched crystal morphologies and absence of typical growth competition gradients in cases of transgranular and mixed type fracturing in small fracture apertures in equally sized and mixed-sized grain structures (Figures 5c, 6a, 8b and 8c).
When the fracture opening velocity is similar to the crystal growth velocity or multiple re-fracturing events happen the vein morphology is impacted by the opening velocity and trajectory (see e.g., Ankit, Selzer, et al., 2015;Ankit, Urai, & Nestler, 2015;Gale et al., 2010;Hilgers et al., 2001;Lander & Laubach, 2015;Prajapati et al., 2018;Urai et al., 1991). Those additional effects are not considered in this work, however in future studies further effects of opening trajectory and multiple (re)fracturing events can be investigated. This will help to better understand the formation of, for example, serrated grain boundaries.
Even tough numerous laboratory studies show calcite growth rates on micro-and macroscale (e.g., Cao et al., 2018;Gratz et al., 1993;Lin & Singer, 2005;Nehrke et al., 2007) very few studies showcase crystal growth in open fractures (e.g., Lee & Morse, 1999) and how it compares to natural fracture opening velocities. Our studies could serve as a pathway for future laboratory experiments (analysis of natural rocks or growth experiments under hydrothermal reservoir conditions) in the direction of precipitation rates of specific facets with the focus on how rates change from rough to faceted growth and how nucleation discontinuities (clay coatings) affect the growth rates.
In the present work, we utilize the same vectors for surface energy (Equation 4) and kinetic anisotropy (Equation 7) and applied the same factor E  (in Equation 7) as in the previous studies of Spruženiece et al. (2020, 2021), which gives good agreement to microscopic observations. When a complete data set from crystal growth experiments for the growth velocities of the crystal facets and the growth velocity drop from rough to euhedral surfaces is available, both the anisotropy vectors and the factor E  (in Equation 7) can be calibrated (as in Prajapati, Abad Gonzalez, et al., 2020;Wendler et al., 2016 for quartz) and quantitative results about the cementation kinetics can be computed. Moreover, for a complete and precise comparison between natural and numerical veins additional data is required, for example: opening vector and shape of fracture surface, pressure, temperature, chemical composition of fluid.
The scaling of the fracture aperture in the simulations by the average host rock diameter m D E enables a direct comparison to natural veins, where the size of the host rock grains differ (e.g., crystalline limestone or mudstone) and helps understanding the conditions, under which the veins have formed. However, in bimodal host rocks, where grains have different sizes (e.g., wackestone, packstone) the scaling by m D E shows different crystal morphologies than in unimodal host rocks (see Figure 9c). For a better comparison of different bimodal host rocks and to unimodal host rocks, where the same mechanisms apply (e.g., fracturing type), the scaling should be extended in future works.
Very few studies (Lander & Laubach, 2015;Lander et al., 2008) have recognized that crystals within a vein can grow with different rates. Our modeling studies are based on prior observational works which inferred that differences in growth rates have to exist (Laubach et al., 2004). As shown by Spruženiece et al. (2020), these growth rate differences can be related to heterogeneities in fracture surface, where transgranularly broken grains grow faster than intergranularly broken grains due to the absence of coatings by secondary minerals which cover the growth surfaces. However, to obtain more quantitative predictions of the models the fracture surface should be investigated more detailed in future works. We simplify the system into fast and slow growing crystals and apply a different mobility for the two different grain classes (factor E ), although in natural fractures multiple grow velocities for IFGs might be present depending on the amount and type of the coating. However, an adapted growth velocity distribution SPÄTH ET AL.

10.1029/2021JB022106
19 of 23 only impacts the grain structure at the fracture rim and does not significantly influence the formation of wider TFGs.
At some point, when the coating is overgrown and crystal facets have formed, the growth velocities between IFGs and TFGs surfaces become similar. For the sake of simplicity, we prescribe a slower growth rate during the whole simulation, even though some IFGs might have formed facets. When this effect is incorporated (time varying E ) similar microstructures form (e.g., wide-blocky veins) to those that form when a higher growth rate ratio E  is prescribed in the early simulation stages, however the crystal morphologies at the rim slightly differ.
The vast majority of our studies were performed in 2D to limit computational costs. In 3D, a higher variation of the crystallographic orientation is possible and the growth of fast-growing crystals is less likely blocked by neighboring crystals as in 2D. However, the 3D simulation shows a similar microstructure as the corresponding 2D simulation (see Figures 9 and 13). This is in agreement with previous studies (Ankit, Selzer, et al., 2015;Ankit, Urai, & Nestler, 2015;Spruženiece et al., 2021;Wendler et al., 2016), which also obtained similar crystal morphologies in their 2D and 3D experiments. We expect that the other investigated parameters also show similar trends in 3D and the findings of our 2D studies can be transferred in 3D as well. Nonetheless, when a porosity and permeability analysis is requested (Kling et al., 2017;Prajapati, Abad Gonzalez, et al., 2020), 3D microstructures are indispensable.
Our results show that vein microstructures can be used as an evidence to infer the conditions of vein formation in nature. In Figure 14 different characteristic crystal morphologies in the simulated veins are compared to natural samples of calcite veins from Somerset coast in southern England and Oman mountains. These crystal morphologies in natural limestone are in accordance with the simulated veins and advocate the capability and accuracy of the phase-field method for modeling vein formation in limestones.
In contrast to the natural veins in mudstone (Figures 14a-14c and 14e) where a diverse crystal morphology exist due to different fracture types (inter-, transgranular, mixed type), the veins in crystalline limestone show mostly blocky crystals (Figures 14c and 2). This indicates that due to the larger grain size either the host rock grains preferentially fractures along weak grain boundaries or that a coating of the host rock grains has a smaller effect and only little fracture surface controlled growth competition occurs ( E   1).

Conclusion and Outlook
The presented simulations demonstrate the capability of the phase-field method to model crystal growth in veins, considering a large number of variables that are expected to influence vein formation in nature. We investigated the effects of fracture surface type, crystal habit, aperture, and grain size distribution in the host rock (analogue to mudstone, wackestone and crystalline limestone in natural limestones) on the vein formation, and created morphology diagrams to represent the results in a systematic way. The simulated vein microstructures show many similarities with natural veins in limestone and overall are in agreement with previous numerical studies (Lander & Laubach, 2015;Nollet et al., 2005;Prajapati et al., 2018).
Depending on the fracture type (intergranular, transgranular or mixed type) and the fracture aperture, different crystal morphologies are observed, such as equidimensional-stretched/ blocky (in small veins), elongated-blocky/ elongated-stretched, and wide blocky.
In future works the findings of this study can be extended into 3D as we showcased with an example. Although no pronounced microstructural differences were observed between 2D and 3D models, the 3D simulations can help to understand the porosity and permeability evolution in veins which is not possible when using the 2D models.
The vein sealing simulations can be combined in future works with mechanical fracture formation models (e.g., Prajapati, Herrmann, et al., 2020) to calculate the morphology of the cracks in porous and fully sealed veins with strong and weak grain boundaries. Furthermore, a coupling of the chemical driving force with SPÄTH ET AL.

Appendix A: Formation of Cross-Shaped Crystals
Here, we demonstrate the formation of cross-shaped crystals which are occasionally found in nature. The fully sealed vein from Section 3.2 with an aperture of 18 m D E and E  = 20 is re-fractured with a high-angle to the first surface ( Figure A1). Similarly, as in Section 3.2 the host rock grains are mostly fractured intergranularly, except two grains in the host rock and the big crystal in the re-fractured vein. Due to the high E  the TFGs expand laterally while overgrowing and terminating the IFGs close to the rim, which results in a cross-shaped crystal.    Figure 3.14 in Arndt (2016)).