Modeling fracture cementation processes in calcite limestone: a phase-field study

Modeling fracture cementation processes in calcite limestone: a phase-field study nishant.prajapati@kit.edu Institute for Applied The present work investigates the influence of crack opening rates on the develop - Materials (IAM-CMS), ment of four important calcite vein morphologies, namely fibrous, elongate-blocky, Karlsruhe Institute of Technology, Building partially open, and euhedral, as a result of bitaxial growth in syntaxial veins using a 30.48, Strasse am Forum 7, multiphase-field model. The continued fracturing that occurs during synkinematic 76131 Karlsruhe, Germany cementation in these veins is simulated using the geometric shift algorithm. The stark Full list of author information is available at the end of the resemblance of the numerically sealed vein microstructures with the natural samples in article terms of structural characteristics as well as remaining pore space signifies a dominant role of crack opening rates in the resulting morphological patterns. Further, simulation results of slow crack opening rates reveal that non-uniform fibers of variable lengths are formed when initial crack aperture is small, due to suppression of growth competi- tion and vice versa. Keywords: Fracture cementation, Multiphase-field model, Crack opening rate Background Veins are ubiquitous geological structures in the earth’s crust and play a pivotal role in deciphering the deformation and fluid-flow histories of the rocks e.g., (Hilgers et  al. 2006; Bons et  al. 2012; McNamara et  al. 2016). In open fractures with circulating for- mation fluids, when temperatures are sufficiently high and geochemical conditions [e.g., pCO , pH, solution ionic composition, saturation state, reacting surface area, presence of inhibiting substances, and solution hydrodynamics (see Plummer et al. 1979; Herman and Lorah 1988; Girard et  al. 1996; Gutjahr et  al. 1996)] favor precipitation of miner- als on the crack surface, syntaxial overgrowth cementation occurs. These fluids serve as the source of vein-filling mineral and are generally supersaturated with respect to the cement-forming solute (Bathurst 1972; Hilgers and Sindern 2005). A wide variety of crack-seal veins are formed due to different combinations of mineral type (e.g., quartz, calcite), growth directions (e.g., syntaxial, antitaxial), number of crack-seal events as well as fracture opening trajectories (e.g., fibrous, blocky, elongate-blocky), as reported in the literature (e.g., Durney and Ramsay 1973; Ramsay 1980; Cox and Etheridge 1983; Lau- bach 1988; Fisher and Brantley 1992; Hilgers and Urai 2002; Laubach et al. 2004a). A deep understanding of crack opening rates, cement growth rates, and vein-forming mechanisms is highly relevant in the analysis of naturally fractured geothermal reser- voirs (e.g., Dezayes et al. 2010; Haffen et al. 2013) as well as hydrocarbon exploration and production (e.g., Gale et  al. 2004; Hilgers et  al. 2006; Lamarche et  al. 2012; Wuestefeld © The Author(s) 2018. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Prajapati et al. Geotherm Energy (2018) 6:7 Page 2 of 15 et  al. 2016). The presence of open cracks enhances flow rates in geothermal projects and hydrocarbon reservoirs in otherwise impermeable formations. The natural frac - ture growth due to repeated fluid pressurization events (i.e., hydrofracturing) signifies a non-continuous and episodic crack opening that might grow rapidly over human time scales. However, recent petrographic investigations of the fractured sandstone in Cre- taceous Travis Peak formation (Becker et  al. 2010) reveal a prolonged fracture growth driven by gas generation, over a duration of 48 m.y. at the estimated crack opening rates of 16–23 μm/m.y. Due to a large number of physical boundary conditions that are present in different geological settings, there exist several mechanisms that could lead to a specific vein tex - ture. For instance, the early work of Taber (1918) proposed that brou fi s veins are formed due to deposition of mineral in the wall rock-vein interface without any fracturing. This growth mechanism was termed as ‘antitaxial’ by Durney and Ramsay (1973), and was later numerically simulated by Urai et  al. (1991) and Hilgers et  al. (2001), and experi- mentally recreated by Means and Li (2001). On the other hand, Mügge (1928) proposed that, for the fractures in their growth stage, if the cement growth rates are faster than the crack opening rates, the growth competi- tion between neighboring crystals is suppressed and brou fi s textures are formed. In gen - eral, the process of cement growth contemporaneous with fracture opening is known as synkinematic cementation (Laubach 1988). Therefore, any cement that has precipitated at the time of fracture opening, whether a crack-seal texture is present or not, is known as a synkinematic deposit. In nature, the sizes of incremental fracture opening can vary from few microns to few millimeters (Fisher and Brantley 1992). Accounts of accurate incremental opening sizes or so-called ‘gaps’ in the aforementioned range have been fur- ther reported in Laubach et al. (2004b) for quartz and in Hooker et al. (2012) for carbon- ate rocks. During synkinematic cementation, a complicated interplay of the continued crack opening and crystal growth occurs. This interplay determines the morphological attributes of the resulting veins based on the crack opening rates (Mügge 1928; Urai et  al. 1991) and was shown to affect quartz-bearing veins (e.g., Laubach 1988; Lander and Laubach 2015), calcite veins (Hilgers et al. 2001), and carbonate veins in dolostones (Gale et  al. 2010). Carbonate vein textures can also be subject to an interplay of syn- and postkinematic cementation (Ukar and Laubach 2016). The analysis of microtextures using e.g., a cathodoluminescence detector mounted on a scanning electron microscope (SEM-CL) might highlight the different cement generations (e.g., Ukar and Laubach 2016). The natural samples of calcite veins derived from the Jurassic limestones of southern England exhibit four prominent morphologies, as depicted in Fig. 1a–d. Different crys - tallographic orientations of vein crystals are reflected by distinct interference colors. The brou fi s veins (in Fig.  1a) are characterized by parallel grain boundaries occurring as thin fibers of different lengths and high length-to-width ratios. Figure  1b, c depicts the vein microstructures with elongate-blocky crystals manifesting crystallographic preferred ori- entations at the final growth stage, which indicates the presence of growth competition during vein formation. Elongate-blocky veins often exhibit partially open morphologies Prajapati et al. Geotherm Energy (2018) 6:7 Page 3 of 15 Fig. 1 Natural samples of calcite veins in Jurassic limestones from southern England exhibiting a fibrous, b elongate-blocky, c partially open, and d euhedral morphologies. The wall rock in fibrous veins is not contained in the sample, but is located towards the base of the image in a. Different interference colors of vein crystals depict different crystallographic orientations. All images are obtained from transmitted microscopy and crossed polarizers [one of the partially open morphologies is also called rind morphologies by Lander and Laubach (2015)]. Here, euhedral crystal faces border some porosity within the partially open fracture (Fig. 1c). Moreover, there exist veins that exhibit euhedral crystals (charac- terized by straight edges and sharp corners when visualized along thin-section) in con- tact with the pore space (black) in open fracture, as illustrated in Fig. 1d. In this work, we utilize a multiphase-field model to delineate the underlying mecha - nisms behind the formation of the aforementioned calcite-vein textures based on the boundary conditions of synkinematic cementation in limestone. The phase-field method is a diffuse interface approach widely used for modeling microstructure evolution pro - cesses such as phase transitions, grain growth, among others (see review articles Boet- tinger et  al. 2002; Moelans et  al. 2008; Nestler and Choudhury 2011; Johannes et  al. 2016). In contrast with the traditional front-tracking methods (see Urai et al. 1991; Bons 2001; Hilgers et  al. 2001; Nollet et  al. 2005), the phase-field approach does not require explicit tracking of the interfaces as the motion of interfaces is inherently captured in the governing equations. This property makes it a computationally efficient and robust approach for the treatment of moving boundary problems such as cement growth in open fractures or porous rocks, even in three dimensions (e.g., Ankit et al. 2013, 2015a, b; Wendler et al. 2016; Prajapati et al. 2017). The present article is organized as follows: In "Methods " section, we describe the phase-field model, the numerical setup and the algorithm for the geometric shift condi - tion. We then present the simulation results of bitaxial calcite cementation in fractures for different crack opening rates and highlight their implications on the resulting vein textures and the histories of fracture porosity in "Results and discussion" section. Fur- ther, the congruence of our numerically simulated vein textures with the natural samples is discussed. In "Conclusion and outlook" section, we conclude the article by recapitulat- ing the important findings and propose directions for further studies. Prajapati et al. Geotherm Energy (2018) 6:7 Page 4 of 15 Methods Multiphase‑field model A thermodynamically consistent multiphase-field model Nestler et  al. (2005) is employed to address the phenomenon of anisotropic crystal growth in calcite veins. The model has been previously applied to address the growth of minerals such as potash alum (Ankit et  al. 2013), quartz (Ankit et  al. 2013, 2015a, b; Wendler et  al. 2016), and calcite (Prajapati et al. 2017). Therefore, for the sake of completeness, we briefly reiterate the model equations. For a detailed discussion and numerical implementation of equa- tions and process-specific model extensions, the interested readers are referred to the above-mentioned articles. We consider a vector of phase-field parameters φ(x, t) =[φ (x, t), . . . , φ (x, t)] , where 1 N φ ∈[0, 1] ∀ α ∈{1, . . . ,N } describes the presence of each calcite grain α at position x and time t. The grain-grain or grain-liquid interface is characterized by a diffuse region of finite width. In the interface region shared between grain α and liquid/other grains, the phase-field parameter φ varies smoothly from 0 outside, to 1 inside the grain. The location and shape of the grain are defined by the isoline with φ = 0.5 . At each spatial point in the system with computational domain  , a summation constraint φ = 1 α=1 is satisfied. The total free energy of the system is formulated as F(φ,∇φ) = εa(φ,∇φ) + w(φ) + f (φ) dΩ , (1) where ε represents a length scale parameter that controls the interfacial-width. εa(φ,∇φ) , w(φ) , and f (φ) are the gradient energy density, the potential energy density, and the bulk free energy density of the system, respectively. The interface region is described by the gradient energy density, given by N ,N 2 2 εa(φ, ∇φ) = ε γ a (φ, ∇φ)|q | , αβ αβ αβ (2) α,β=1 (α<β) where γ denotes the interfacial energy density between α and β grains, and q αβ αβ (= φ ∇φ − φ ∇φ ) represents the generalized gradient vector pointing in the outward α β β α normal direction of the α–β interface. The choice of function a (φ,∇φ) determines the αβ type of interface energy anisotropy ( a = 1 for the isotropic case). Although, calcite αβ exhibits a wide variety of growth habits including acute to obtuse rhombohedral, sca- lenohedral, prismatic, among several others, we limit the present discussion to a com- mon rhombohedral calcite geometry (see Fig. 2a). In the present study, we assume that the crystal facets are developed due to strong faceted anisotropy in the interface energy. Therefore, a piece-wise function of the form αβ a (φ, ∇φ) = max · η αβ (3) k,αβ 1≤k≤n |q | αβ αβ is utilized for incorporating such faceted anisotropy in the interface energy formula- tion. Here, {η | k = 1, . . . , n } denotes the position vectors corresponding to the n k,αβ αβ αβ Prajapati et al. Geotherm Energy (2018) 6:7 Page 5 of 15 Right−handed Symmetric Left−handed Diffuse 1 interface liquid 90° −90° Top views φ grain grain φ =0.5 grain t=0 t=1505 a Initial spherical grain Simulated shape Rhombohedral calcite geometry Fig. 2 a Symmetric and asymmetric projections of the rhombohedral calcite geometry. b Validation of interface energy anisotropy formulation corresponding to the asymmetric left-handed 2D projection. A spherical grain evolves to its equilibrium shape while retaining a constant volume. The final crystal shares a diffuse interface of finite width with the surrounding liquid. The shape and size of the crystal are defined by the isoline corresponding to φ = 0.5 grain Table 1 Values of phase-field model parameters used in simulations Model parameters Value Grid size x 1 Time-step width t 0.07 Length scale parameter ε 6.5 Interfacial energy density γ , i.e., γ and γ 1 αβ g,g g,l Higher order parameter γ 10 αβδ Kinetic coefficient for grain–grain interface τ 100 g,g Kinetic coefficient for grain–liquid interface τ 1 g,l vertices of the Wulff shape of crystal α dispersed in bulk phase β . The evolution equa - tions of the phase-fields read ∂φ ∂a(φ, ∇φ) ∂a(φ, ∇φ) 1 ∂w(φ) ∂f (φ) τε =ε ∇· − − − − (4) ∂t ∂∇φ ∂φ ε ∂φ ∂φ α α α α for α = 1, . . . ,N . Here, τ is the interfacial kinetic coefficient and  is the Lagrange parameter to ensure the summation constraint φ = 1 . The set of partial differential equations in Eq. α=1 (4) is solved using forward Euler scheme for time derivative and second-order accurate central difference scheme for the spatial derivatives. The numerical solution method is efficiently implemented in a parallel simulation framework for crystal evolution, PACE3D (see Hötzer et al. 2018), which is written in C language. The values of model parameters used in the present work are listed in Table 1. Modeling calcite cements The present work comprises a two-dimensional treatment of the complex process of synkinematic cementation with continuous fracture growth. Therefore, we obtain the following 2D projections: right-handed, symmetric and left-handed (see Fig.  2a), that were also highlighted in our previous article pertaining to three-dimensional mode- ling of cement overgrowth in porous limestones (Prajapati et al. 2017). During bitaxial Prajapati et al. Geotherm Energy (2018) 6:7 Page 6 of 15 cement growth in fractures, the effect of asymmetricity can be significant (Ankit et  al. 2015b). Hence, we chose the left-handed projection of the calcite, in order to account for the influence of asymmetricity. The chosen position vectors for the vertices of crystal shape are as follows: η = [0, ±0.75] η = [1, 0.25]= −η 1,2 3 4 resulting in an aspect ratio of 1.33 of the 2D crystal. Figure  2b depicts the simulated left-handed 2D equilibrium shape of calcite from an initial spherical grain embedded in liquid phase, using the volume preservation technique (Nestler et al. 2008). For sim- ulating bitaxial cement growth in progressive fracture-opening, the numerical setup is described in "Numerical setup" section. Numerical setup The initial grain structure in a 2D computational domain of dimensions 1000x × 1000x is generated using a voronoi diagram based on a random set of points, see Fig. 3a. The gray region below the lower crack surface is rendered stationary and is termed as barrier phase. In the very first step, a fracture of finite width is intro - duced resulting in two crack surfaces with randomly oriented grains as shown in Fig. 3b. The fracture is completely filled with liquid phase (in yellow), which is assumed to be supersaturated with respect to calcite. This assumption is realized by applying a constant driving force ( f ) for the anisotropic crystal growth. In order to repress the influence of the opposing curvature-dependent forces on the crystal growth (see Plapp (2012) for fundamentals of phase-field), a sufficiently high value of f (=− 0.15) is chosen for all the simulations. The defined boundary conditions for the numerical studies are shown in Fig.  3c. The geometric shift boundary condition applied at the lower crack surface is described in "Geometric shift condition" section. The definition of crystal orientations is depicted in Fig. 3d. Geometric shift condition In order to simulate an opening crack, the previous unitaxial crack opening algorithm utilized by Ankit et  al. (2013) was amended for bitaxial growth in syntaxial veins. The present so-called geometric shift algorithm relocates the entire set of grains on the lower surface by a given vertical displacement after a prescribed number of time steps and fills a b c Isolate d 1000 Δ x Liquid Voronoi crack Geometric shift Periodic Periodic 50 Δ x surfaces condition with grains Initial crack introduced orientation Initial Barrier crack width −90 Fig. 3 a Dimensions of the computational domain consisting of a layer of polycrystalline grains (in RGB colors) generated using voronoi diagram and the stationary barrier phase (in gray). b An initial crack filled with liquid (in yellow) is introduced. c The applied boundary condition at the edges of computational domain. d Definition of crystallographic orientation of calcite grains 1000 Δ x Prajapati et al. Geotherm Energy (2018) 6:7 Page 7 of 15 Table 2 Values of  the  prescribed initial crack width and  the  crack opening rates for  three different cases Cases Initial crack width Crack opening rate (COR) Slow opening rate 100x �x/1000�t Moderate opening rate 100x �x/200�t Fast opening rate 100x �x/100�t the new empty cells with liquid phase. This essentially implies that, in the scenario of crystal-bridging or vein closure, the crack always opens along the boundaries between the grains that lie on mutually opposite crack surfaces. Moreover, the phase-field evolu - tion is not computed in the barrier phase, but the barrier is constantly consumed due to the shifting process. Results and discussion Eec ff t of crack opening rate: implications for vein morphologies By employing the generated computational domain with randomly oriented calcite seeds on the crack surfaces (in Fig.  3a), and the geometric shift condition for crack opening (described in "Geometric shift condition" section), simulations were performed for equal initial crack width and different values of crack opening rates, see Table  2 for param- eters. Each simulation was performed until the barrier phase was completely consumed due to the moving lower crack surface along with the attached grains. In particular, three cases of slow, moderate, and fast opening rates were considered as follows: • Slow opening rate: For the applied driving force ( f =− 0.15 ), the prescribed crack opening rate (in Table 2) is smaller than the minimum potential growth. For the cho- sen crack width of 100x (Fig. 4a), the magnitude of initial crack aperture is one order higher than the initial mean grain size of aggregates. This leads to a growth compe - tition in the early stages as illustrated at t = 1190 in Fig.  4b. As soon as the crystals grow and fill the initial crack (at t = 3850 in Fig.  4b), further growth of crystals is almost completely dependent on the opening increments. From this stage onwards, the growth competition arising due to random grain orientations is suppressed. This implies that almost all grains grow equally into the fracture spaces introduced dur- ing small crack opening increments. Different stages of synkinematic cementation are shown in Fig. 4a, b. The resulting vein consists of brou fi s crystals of nearly equal width and parallel to the crack opening displacements (at t = 59605 in Fig. 4b, e). • Moderate opening rate: For the same driving force of grain growth as the previous case, the prescribed crack opening rate (in Table  2) lies between the maximum and minimum potential growth rate. Different stages of synkinematic cementation are shown in Fig. 4c. In the initial stages, due to the availability of sufficient pore space, the growth competition between the neighboring crystals of the same crack surface is dominant. It is noteworthy that since the aspect ratio of the chosen 2D equilibrium shape of calcite is close to unity, the orientation selection is more or less random. Therefore, the fastest growing crystals are decided on the basis of neighboring grain Prajapati et al. Geotherm Energy (2018) 6:7 Page 8 of 15 Early growth competition Suppressed growth due to high initial competition after crack−width first complete sealing Slow crack opening rate Fibers t=1190t=3850 t=59605 e crack opening direction Bridging crystals Moderate crack opening Partially−open rate vein Elongate−blocky t=0.7 t=2660 t=4935 t=11935 f a c Fast crack opening rate Open vein at all stages Euhedral t=1050 t=2450 t=5985 d g Fig. 4 a Computational domain with an initial crack width of 100x . Simulated bitaxial crystal growth for different values of crack opening rates: b Slow crack opening rate �x/1000�t, c moderate crack opening rate �x/200�t, and d fast crack opening rate �x/100�t . The progression is shown at representative simulation time steps. The selected crystals in the resulting vein morphologies are illustrated. e Fibers with almost equal width. f Elongate-blocky crystals with increasing width. g Crystals forming euhedral flat facets due to availability of sufficient pore space for crystal growth interactions. As the cementation progresses, the crystals begin to form bridges (at t = 2660 in Fig.  4c), so that the further growth is constrained by the crystals grow- ing from the opposite crack surface resulting in partially open veins (at t = 4935 in Fig.  4c). At later stages, the pore space is entirely filled with bitaxial crystal growth leading to elongate-blocky vein textures with increasing width along the length of the crystals (at t = 11935 in Fig. 4c, f ). • Fast opening rates: In this case, the crack opening rates are higher than the maxi- mum potential growth rate. Under such conditions, the growth competition between the neighboring crystals is present and the growth is not constrained by the crys- tals growing from the opposite fracture surface. The veins are open at all stages of time and the resulting crystals are euhedral in nature with flat facets and sharp cor - ners, as depicted at t = 5985 in Fig. 4d, g. The texture resembles rind morphologies described in Lander and Laubach (2015). In general, the aforementioned simulated results for calcite cement growth in fractures, depending upon fracture opening rate and crystal orientations are in good agreement with natural and experimental results, presented previously for quartz cements in fractures (Lander and Laubach 2015). Eec ff t of crack opening rate: implications on porosity evolution The phase-field simulations (in "Effect of crack opening rate: implications for vein morphologies" section) provide valuable insights into the processes controlling the Prajapati et al. Geotherm Energy (2018) 6:7 Page 9 of 15 formation of fibrous, partially open, elongate-blocky, and euhedral veins based on the boundary conditions of synkinematic cementation. Apart from recreating the mech- anisms behind various vein morphologies, the numerical simulations can provide information on histories of fracture porosity as a result of synkinematic cementation for different crack opening rates, see Fig.  5. The porosity-time behavior results from the interplay of crack opening rate and the anisotropic growth of crystals. For the present study, constant crack opening rates are applied for all cases, thereby implying a linear porosity-time behavior in the absence of anisotropic crystal growth (plotted in same colored dashed lines as the fracture sealing curves in Fig. 5). The anisotropic growth of crystals adds the non-linearity in the time evolution of porosity for the three crack opening rates (slow, moderate and fast). We analyze each case one by one as follows: • Slow opening rate (crack opening rate = �x/1000�t ): The porosity-time curve monotonically decreases with respect to time. Due to the early growth competition (described in "Effect of crack opening rate: implications for vein morphologies " sec- tion), the initial decline in the porosity is slow but accelerates as the cementation proceeds. After the first complete sealing (somewhere close to t = 3850 as shown in Fig. 4b), the crack reopening occurs along the interfaces between the crystals grow- ing from opposite sides of the vein. The subsequent porosity evolution is completely dependent on the crack opening increment. The temporal evolution of porosity enters a periodic oscillatory regime due to the occurrence of alternate events of crack opening and vein filling, as illustrated in the zoomed inset picture of Fig. 5. • Moderate opening rate (crack opening rate = �x/200�t ): The porosity initially increases with time owing to the early growth competition (that suppresses crystal growth) and a higher opening rate. The growth competition declines due to occlu - sion of the growth of unfavorably oriented crystals, resulting in an increase in the net crystal growth rates. These increased growth rates lead to a maximum (at t = 1225 ) followed by a monotonic decrease in the porosity-time curve. Due to the higher crack opening rate than the former case, the first complete sealing takes place at later stage in the present case. 0.8 Δx/ 1000Δt Δx/ 200Δt 0.7 Δx/ 100Δt Δx/ 50 Δt t=5985 0.6 t=3010 0.5 0.4 t=59605 0.3 t=0.7 0.2 Periodic 0.131 0.1 t=11935 0 2000 4000 6000 8000 10000 12000 59605 1225 time Fig. 5 Temporal evolution of fracture porosity for different crack opening rates represented by different colored solid lines. The same colored dashed lines depict the porosity-time behavior in the absence of cement growth, signifying constant crack opening rates. The inset picture illustrates the initial and final vein morphologies at the end of simulation porosity Prajapati et al. Geotherm Energy (2018) 6:7 Page 10 of 15 • Fast opening rate (crack opening rates = �x/100�t and �x/50�t ): The evolution of porosity monotonically increases with time. For even larger opening rates (i.e., �x/50�t ), the porosity-time curve gets steeper. It is noteworthy that the porosity- time curve has an upper bound that is prescribed by the crack opening rates (dashed lines in Fig. 5). Therefore, although the porosity increases monotonically, the poros - ity-time curve is concave with respect to time. Eec ff t of initial crack aperture From the previous set of simulations (in "Effect of crack opening rate: Implications for vein morphologies" section), it is evident that the growth competition of different grain orientations is suppressed by the slow crack opening rates leading to fibrous morphol - ogies. However, the early growth competition that arises due to the large initial crack aperture (one order higher than the grain size) results in the consumption of less favora- bly oriented grains in the initial stages. In this section, we investigate the influence of initial crack width on the resulting fibrous morphologies. In particular, we consider two cases of small and large initial cracks. The parameters used for the present study are listed in Table 3. – Small initial crack: In this case, the magnitude of initial crack width and grain size are of same order. Under such conditions, the growth competition is suppressed from the very beginning, leading to retardation in the occlusion of unfavorably oriented grains by the favorable ones. Therefore, more initial grains are preserved as shown at t = 8400 in Fig. 6a. As the cementation proceeds, the unfavorably oriented crystals are gradually occluded after the initial stages, resulting in arrested fibers of reduced length at intermediate time step ( t = 31150 in Fig.  6a). The resulting veins (at t = 66220 in Fig. 6a) possess non-uniform fibrous textures consisting of fibers of different lengths. – Large initial crack: The magnitude of crack width is one order higher than that of grain size. Therefore, the intense initial growth competition (at t = 1890 in Fig.  6b) leads to occlusion of almost all the unfavorably oriented grains in the early stages. The fib - ers formed as a result of growth of favorably oriented crystals during the intermediate stages (at t = 24535 in Fig.  6b) exhibit very little growth competition. Therefore, the resulting veins (at t = 59605 in Fig. 6b) have uniform fibrous morphologies composed of almost all fully grown non-occluded fibers. In the following sections, we refer to these fibers as the ‘survivor’ fibers that have continued to grow from the crack surface to the center of the vein, without being occluded by the neighboring crystals. We plot the temporal evolution of porosity along with the initial and final vein mor - phologies for the above cases in Fig.  7a. For a small initial crack width, the opening Table 3 Values of the prescribed initial crack width and the crack opening rate Cases Initial crack width Crack opening rate Small initial crack 5x �x/1000�t Large initial crack 100x �x/1000�t Prajapati et al. Geotherm Energy (2018) 6:7 Page 11 of 15 t=66220 Survivor fibers Small crack t=31150 t=8400 t=0.7 ICW: 5 Δ x Negligible growth More intial grain competition orientations are preserved Arrested fibers VW: 177 Δ x VW: 500 Δ x VW: 1000 Δ x Almost all survivor fibers Large crack b Initial growth Fewer initial grain orientations are competition preserved ICW: 100 Δ x t=0.7 t=1890 t=24535 ICW: Initial crack width Increasing vein width vein−width (VW) t=59605 Fig. 6 Simulated bitaxial crystal growth for a initial crack width = 5x (small crack) and b initial crack width = 100x (large crack). Progression is shown at representative values of the evolving vein-width ( VW ). For small cracks, the resulting veins (at t = 66220 in a) consist of arrested as well as survivor fibers of different lengths, while large cracks lead to uniform fibrous vein textures (at t = 59605 in b) 0.12 70 Small initial crack Small initial crack crack opening rate width (=5 Δx) Δx) width (=5 Δx/ 1000Δt 0.1 Large initial crack Large initial crack width (=100 Δx) width (=100 Δx) 0.08 t=0.7 t=59605 t=66220 0.06 0.04 Periodic 0.02 0 20 200 300 400 500 600 700 800 900 1000 0 2000 4000 60000 66220 time vein width [ x] Δ a t=0.7 b Fig. 7 a Temporal evolution of fracture porosity for different initial crack widths. The inset pictures illustrate the initial and final vein morphologies at the end of simulation. b Plot of the number of survivor grains with the evolving vein-width for different values of initial crack width increments are of same order, forcing the cement growth to be entirely dependent on the crack opening. Therefore, the porosity evolution enters the periodic regime at an ear - lier stage than for the larger crack opening. As highlighted in Fig. 6, the simulated over- growths for the two cases are compared at same vein-width (VW). In order to analyze the influence of initial growth competition, we plot the number of survivor grains with the evolving vein-width, see Fig. 7b. For a vein-width of 200x , the number of survivor grains is lower for the large ini- tial crack due to the suppression of growth of unfavorably oriented grains in the very beginning. The remaining grains grow with little competition; thereby, the decline in the number of survivor grains is lower as compared to the case of small initial crack. For the small initial crack, the unfavorably oriented grains are not occluded initially and are reduced gradually over the entire process of cementation. It is noteworthy that the resulting number of survivor fibers is similar in both cases. Present simulations provide valuable insights into the growth dynamics of arrested and survivor fibers. We infer that natural samples of non-uniform fibrous veins exhibiting arrested and survivor fibers porosity number of survivor grains Prajapati et al. Geotherm Energy (2018) 6:7 Page 12 of 15 result due to slow opening rates and small crack apertures at all stages of synkinematic cementation. On the other hand, uniformly long fibers in the fibrous veins indicate the presence of large crack apertures at the initial stages of synkinematic cementation, lead- ing to the occlusion of unfavorably oriented crystals in the very beginning. Comparison with natural samples Natural samples of calcite veins showcase a wide variety of vein cement morpholo- gies, which can be linked to the simulation results. All shown samples were observed in Jurassic limestones from southern England (Blue Lias at Kilve and Lyme Regis, hand specimens were collected from loose boulders at the beaches). Euhedral vein-form- ing crystals, characterized by straight edges in contact with open pore space (black in Fig.  8a), are in good agreement with simulated examples, highlighting growth compe- tition at the contact with the host rock and increasing crystal width towards the pore (Fig. 8a). Pore morphologies and occasional irregular domain contacts between bitaxially grown calcite vein cements in partially open veins show a good match between natural and simulated examples (Fig. 8b). This includes the possibility of euhedral crystal facets in small open pores. The elongate-blocky morphology of moderate crack opening rates, resulting in a mostly sealed vein is characterized by irregular contacts between domains grown from opposing sides of the crack (Fig. 8c). This also includes growth competition at the host rock interface of the fracture. The case of slow crack opening rate and small initial crack aperture resulting in non-uniform fibrous calcite cement of varying length matches well with habits observed in natural brou fi s veins (Fig. 8d). Conclusion and outlook While the recent numerical works (e.g., Ankit et al. 2015a, b; Wendler et al. 2016) advo- cate the potential of the phase-field approach in addressing quartz mineral growth in sandstones, the present study in addition to our previous work (Prajapati et  al. 2017) serves as a basis for new calcite cementation models for carbonate rocks based on sound thermodynamic principles. The present study complements the works of Bons (2001) and Hilgers et  al. (2001) (who used front-tracking technique) and adds to them. While the vein morphologies arising due to the interplay of different crack opening rates and Fast crack−opening rate Moderate crack−opening rate Slow crack−opening rate a Euhedral veins Partially−openaElongate−blocky d Fibrous b c suppressed grains Euhedral crystal facets in small open spaces suppressed grains due to growth Euhedral crystals Vein fibers of different lengths competition at host rock interface Fig. 8 Comparison of simulated vein morphologies with natural samples. The sketches of natural samples (in gray scale) are drawn from the natural samples of calcite veins in Jurassic limestones from southern England, shown in Fig. 1. a Euhedral, b partially open, c elongate-blocky and d fibrous veins. The wall rock (gray) in the sketch of fibrous vein is not contained in the sample, but it is towards the base of the image. The simulated veins exhibit strong resemblance with the natural samples suggesting a pronounced role of crack opening rates and initial crack width in the formation of different vein textures Prajapati et al. Geotherm Energy (2018) 6:7 Page 13 of 15 cement growth rates are faithfully reproduced by the phase-field simulations, the present study additionally highlights the formation mechanisms of uniform and non-uniform fibrous veins based on the initial crack apertures. The convincing resemblance between the numerically sealed vein morphologies with natural samples further demonstrates the capabilities of the employed multiphase-field model in elucidating the growth dynam - ics of four prominent calcite-vein textures (fibrous, elongate-blocky, partially open, and euhedral). Furthermore, the model in conjunction with innovative pre- and post-pro- cessing tools developed for the analysis and visualization of numerical data allows to speculate the reduction of fracture porosity over time. Although the present study is limited to 2D, the implementation of the model in 3D is fairly straight forward, albeit with increased computational resources. This provides a distinct edge over the conventional front-tracking methods that are restricted to 2D owing to numerical issues that arise while addressing the process in three dimensions. Although bitaxial crack-seal mechanism is by far the most common vein-forming pro- cess, nevertheless, very few numerical studies have hitherto been conducted, due to the aforementioned reasons. The extension of the present work to three dimensions remains an interesting prospect for future studies that enables the analyses of the fracture net- works and flow behavior in the geothermal reservoir rocks. The present work assumes a constant crack opening and models fracturing as a continuous process. However, the implemented geometric shift condition can directly be utilized to simulate episodic reopening of fractures (orthogonal to the crack surface) in response to repeated pres- surization events, i.e., hydrofracturing. The present work paves the way for the ongoing studies that explore the development of syntaxial textures due to complex fracture open- ing histories and paths such as crack opening obliquely to the vein boundaries or along curved trajectories. Moreover, building upon the previous works of (Ankit et al. 2013, b), the model along with appropriate extensions can be utilized to further explore the controls on the formation of antitaxial and ataxial veins. The present work is an initial treatment of a highly coupled thermo-hydro-mechanical-chemical problem of vein for- mation as a result of cementation, wherein faceted crystal growth mechanisms based on the boundary conditions of synkinematic cementation are simulated. The present phase- field model needs to be extended in order to account for the mechanical forces for crack opening and chemical driving forces for crystal growth that facilitates the differential growth rates on facetted and unfacetted surfaces as suggested by Ukar and Laubach (2016). Provided the availability of physical parameters (e.g., crystal orientations, super- saturation of solute, absolute rates of calcite accumulation) through benchmark cement growth experiments, the model can be calibrated to precisely mimic the experimentally synthesized veins and further explore the controls on fabric development in natural veins. Authors’ contributions NP conducted the numerical studies under the supervision of MS in addition toconceptual discussions with BN, BB, and CH. MS and NP implemented the geometric shift algorithm for the present study. BB and CH provided the micropho- tographs of the natural samples. NP drafted the manuscript. All authors editedthe manuscript. All authors read and approved the final manuscript. Author details Institute for Applied Materials (IAM-CMS), Karlsruhe Institute of Technology, Building 30.48, Strasse am Forum 7, 76131 Karlsruhe, Germany. Institute of Digital Materials Science (IDM), Karlsruhe University of Applied Sciences, Molt- kestr. 30, 76133 Karlsruhe, Germany. Institute of Applied Geosciences (AGW-SGT ), Karlsruhe Institute of Technology, Adenauerring 20a, 76131 Karlsruhe, Germany. Prajapati et al. Geotherm Energy (2018) 6:7 Page 14 of 15 Acknowledgements The authors would like to thank the Helmholtz Association for the financial support through “KIT-Integration initiative in geothermal systems” and within the program “RE-renewable energies.” Further, we acknowledge publication support by Helmholtz Centre for Environmental Research—UFZ; Helmholtz Centre Potsdam—GFZ German Research Centre for Geosciences; and Karlsruhe Institute of Technology for the article processing charges. Competing interests The authors declare that they have no competing interests. Availability of data and materials All relevant data and material are presented in the paper. Consent for publication Not applicable. Ethics approval and consent to participate Not applicable. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Received: 15 March 2018 Accepted: 26 May 2018 References Ankit K, Nestler B, Selzer M, Reichardt M. Phase-field study of grain boundary tracking behavior in crack-seal microstruc- tures. Contrib Miner Petrol. 2013;166(6):1709–23. Ankit K, Selzer M, Hilgers C, Nestler B. Phase-field modeling of fracture cementation processes in 3-D. J Petrol Sci Res. 2015;4(2):79–96. Ankit K, Urai JL, Nestler B. Microstructural evolution in bitaxial crack-seal veins: a phase-field study. J Geophys Res. 2015;120(5):3096–118. Bathurst RG. Carbonate sediments and their diagenesis, vol. 12. Amsterdam: Elsevier; 1972. Becker SP, Eichhubl P, Laubach SE, Reed RM, Lander RH, Bodnar RJ. A 48 my history of fracture opening, temperature, and fluid pressure: Cretaceous Travis Peak Formation. East Texas basin. Bulletin. 2010;122(7–8):1081–93. Boettinger WJ, Warren JA, Beckermann C, Karma A. Phase-field simulation of solidification. Ann Rev Mater Res. 2002;32(1):163–94. Bons PD. Development of crystal morphology during unitaxial growth in a progressively widening vein: I. The numerical model. J Struct Geol. 2001;23(6–7):865–72. Bons PD, Elburg MA, Gomez-Rivas E. A review of the formation of tectonic veins and their microstructures. J Struct Geol. 2012;43:33–62. Cox SF, Etheridge MA. Crack-seal fibre growth mechanisms and their significance in the development of oriented layer silicate microstructures. Tectonophysics. 1983;92(1–3):147–70. Dezayes C, Genter A, Valley B. Structure of the low permeable naturally fractured geothermal reservoir at Soultz. Comptes Rendus Geosci. 2010;342(7–8):517–30. Durney D, Ramsay JG. Incremental strains measured by syntectonic crystal growths. In: De Jong KA, Scholten R, editors. Gravity and tectonics. New York: Wiley; 1973. p. 67–96. Fisher DM, Brantley SL. Models of quartz overgrowth and vein formation: deformation and episodic fluid flow in an ancient subduction zone. J Geophys Res. 1992;97(B13):20043–61. Gale JF, Laubach SE, Marrett RA, Olson JE, Holder J, Reed RM. Predicting and characterizing fractures in dolostone reser- voirs: using the link between diagenesis and fracturing. Geol Soc Lond. 2004;235(1):177–92. Gale JF, Lander RH, Reed RM, Laubach SE. Modeling fracture porosity evolution in dolostone. J Struct Geol. 2010;32(9):1201–11. Girard JP, Sanjuan B, Czernichowski-Lauriol I, Fouillac C. Diagenesis of the Oseberg Sandstone Reservoir (North Sea): an example of integration of core, formation fluid and geochemical modelling studies. AAPG Bulletin. 5(CONF-960527). Gutjahr A, Dabringhaus H, Lacmann R. Studies of the growth and dissolution kinetics of the C aCO polymorphs calcite and aragonite II. The influence of divalent cation additives on the growth and dissolution rates. J Cryst Growth. 1996;158(3):310–5. Haffen S, Géraud Y, Diraison M, Dezayes C. Fluid-flow zones in a geothermal sandstone reservoir: localization from ther - mal conductivity and temperature logs, borehole EPS1(Soultz-sous-Forêts, France). Geothermics. 2013;46:32–41. Herman JS, Lorah MM. Calcite precipitation rates in the field: measurement and prediction for a travertine-depositing stream. Geochimica et Cosmochimica Acta. 1988;52(10):2347–55. Hilgers C, Koehn D, Bons PD, Urai JL. Development of crystal morphology during unitaxial growth in a progressively wid- ening vein: II. Numerical simulations of the evolution of antitaxial fibrous veins. J Struct Geol. 2001;23(6–7):873–85. Hilgers C, Sindern S. Textural and isotopic evidence on the fluid source and transport mechanism of antitaxial fibrous microstructures from the Alps and the Appalachians. Geofluids. 2005;5(4):239–50. Hilgers C, Urai JL. Microstructural observations on natural syntectonic fibrous veins: implications for the growth process. Tectonophysics. 2002;352(3–4):257–74. Prajapati et al. Geotherm Energy (2018) 6:7 Page 15 of 15 Hilgers C, Nollet S, Schonherr J, Urai JL. Paleo-overpressure formation and dissipation in reservoir rocks. Oil Gas Eur Mag. 2006;32(2):68–73. Hilgers C, Kirschner DL, Breton JP, Urai JL. Fracture sealing and fluid overpressures in limestones of the Jabal Akhdar dome. Oman mountains. Geofluids. 2006;6(2):168–84. Hooker JN, Gomez LA, Laubach SE, Gale JFW, Marrett R. Eec ff ts of diagenesis (cement precipitation) during fracture open- ing on fracture aperture-size scaling in carbonate rocks. Geol Soc Lond. 2012;370(1):187–206. Hötzer J, Reiter A, Hierl H, Steinmetz P, Selzer M, Nestler B. The parallel multi-physics phase-field framework Pace3D. J Comput Sci. 2018;26:1–12. Johannes H, Michael K, Philipp S, Britta N. Applications of the phase-field method for the solidification of microstructures in multi-component systems. J Indian Inst Sci. 2016;96(3):235–56. Lamarche J, Lavenu AP, Gauthier BD, Guglielmi Y, Jayet O. Relationships between fracture patterns, geodynamics and mechanical stratigraphy in Carbonates (South-East Basin, France). Tectonophysics. 2012;581:231–45. Lander RH, Laubach SE. Insights into rates of fracture growth and sealing from a model for quartz cementation in frac- tured sandstones. GSA Bull. 2015;127(3–4):516–38. Laubach SE. Subsurface fractures and their relationship to stress history in East Texas Basin sandstone. Tectonophysics. 1988;156(1–2):37–49. Laubach SE, Reed RM, Olson JE, Lander RH, Bonnell LM. Coevolution of crack-seal texture and fracture porosity in sedi- mentary rocks: cathodoluminescence observations of regional fractures. J Struct Geol. 2004;26(5):967–82. Laubach SE, Lander RH, Bonnell LM, Olson JE, Reed RM. Opening histories of fractures in sandstone. Geol Soc Lond. 2004;231(1):1–9. McNamara DD, Lister A, Prior DJ. Calcite sealing in a fractured geothermal reservoir: insights from combined EBSD and chemistry mapping. J Volcanol Geotherm Res. 2016;323:38–52. Means WD, Li T. A laboratory simulation of fibrous veins: some first observations. J Struct Geol. 2001;23(6–7):857–63. Moelans N, Blanpain B, Wollants P. An introduction to phase-field modeling of microstructure evolution. Calphad. 2008;32(2):268–94. Mügge O. Ueber die Entstehung faseriger Minerale und ihrer Aggregationsformen. Neues Jahrbuch für Mineralogie, Geologie und Paläontologie A. 1928;58:303–438. Nestler B, Choudhury A. Phase-field modeling of multi-component systems. Curr Opin Solid State Mater Sci. 2011;15(3):93–105. Nestler B, Garcke H, Stinner B. Multicomponent alloy solidification: phase-field modeling and simulations. Phys Rev E. 2005;71(4):041609. Nestler B, Wendler F, Selzer M, Stinner B, Garcke H. Phase-field model for multiphase systems with preserved volume frac- tions. Phys Rev E. 2008;78(1):011604. Nollet S, Urai JL, Bons PD, Hilgers C. Numerical simulations of polycrystal growth in veins. J Struct Geol. 2005;27(2):217–30. Plapp Mathis. Phase-field models. Multiphase microfluidics: the diffuse interface model. Vienna: Springer; 2012. p. 129–75. Plummer LN, Parkhurst DL, Wigley TML. Critical review of the kinetics of calcite dissolution and precipitation. In: Jenne E. Chemical modeling in aqueous systems; ACS Symposium Series. 1979, p. 537-73. Prajapati N, Selzer M, Nestler B. Computational modeling of calcite cementation in saline limestone aquifers: a phase- field study. Geotherm Energy. 2017;5(1):15. Ramsay JG. The crack-seal mechanism of rock deformation. Nature. 1980;284:135–9. Taber S. The origin of veinlets in the Silurian and Devonian strata of central New York. J Geol. 1918;26(1):56–73. Ukar E, Laubach SE. Syn-and postkinematic cement textures in fractured carbonate rocks: insights from advanced cathodoluminescence imaging. Tectonophysics. 2016;690:190–205. Urai JL, Williams PF, Van Roermund HLM. Kinematics of crystal growth in syntectonic fibrous veins. J Struct Geol. 1991;13(7):823–36. Wendler F, Okamoto A, Blum P. Phase-field modeling of epitaxial growth of polycrystalline quartz veins in hydrothermal experiments. Geofluids. 2016;16(2):211–30. Wuestefeld P, De Medeiros M, Koehrer B, Sibbing D, Kobbelt L, Hilgers C. Automated workflow to derive LIDAR fracture statistics for the DFN modelling of a tight gas sandstone reservoir analog. In: 78th EAGE Conference and Exhibition 2016; 2016. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geothermal Energy Springer Journals

Modeling fracture cementation processes in calcite limestone: a phase-field study

Free
15 pages
Loading next page...
 
/lp/springer_journal/modeling-fracture-cementation-processes-in-calcite-limestone-a-phase-szIcjlxVzR
Publisher
Springer Berlin Heidelberg
Copyright
Copyright © 2018 by The Author(s)
Subject
Earth Sciences; Geotechnical Engineering & Applied Earth Sciences; Renewable and Green Energy; Geoecology/Natural Processes
eISSN
2195-9706
D.O.I.
10.1186/s40517-018-0093-4
Publisher site
See Article on Publisher Site

Abstract

nishant.prajapati@kit.edu Institute for Applied The present work investigates the influence of crack opening rates on the develop - Materials (IAM-CMS), ment of four important calcite vein morphologies, namely fibrous, elongate-blocky, Karlsruhe Institute of Technology, Building partially open, and euhedral, as a result of bitaxial growth in syntaxial veins using a 30.48, Strasse am Forum 7, multiphase-field model. The continued fracturing that occurs during synkinematic 76131 Karlsruhe, Germany cementation in these veins is simulated using the geometric shift algorithm. The stark Full list of author information is available at the end of the resemblance of the numerically sealed vein microstructures with the natural samples in article terms of structural characteristics as well as remaining pore space signifies a dominant role of crack opening rates in the resulting morphological patterns. Further, simulation results of slow crack opening rates reveal that non-uniform fibers of variable lengths are formed when initial crack aperture is small, due to suppression of growth competi- tion and vice versa. Keywords: Fracture cementation, Multiphase-field model, Crack opening rate Background Veins are ubiquitous geological structures in the earth’s crust and play a pivotal role in deciphering the deformation and fluid-flow histories of the rocks e.g., (Hilgers et  al. 2006; Bons et  al. 2012; McNamara et  al. 2016). In open fractures with circulating for- mation fluids, when temperatures are sufficiently high and geochemical conditions [e.g., pCO , pH, solution ionic composition, saturation state, reacting surface area, presence of inhibiting substances, and solution hydrodynamics (see Plummer et al. 1979; Herman and Lorah 1988; Girard et  al. 1996; Gutjahr et  al. 1996)] favor precipitation of miner- als on the crack surface, syntaxial overgrowth cementation occurs. These fluids serve as the source of vein-filling mineral and are generally supersaturated with respect to the cement-forming solute (Bathurst 1972; Hilgers and Sindern 2005). A wide variety of crack-seal veins are formed due to different combinations of mineral type (e.g., quartz, calcite), growth directions (e.g., syntaxial, antitaxial), number of crack-seal events as well as fracture opening trajectories (e.g., fibrous, blocky, elongate-blocky), as reported in the literature (e.g., Durney and Ramsay 1973; Ramsay 1980; Cox and Etheridge 1983; Lau- bach 1988; Fisher and Brantley 1992; Hilgers and Urai 2002; Laubach et al. 2004a). A deep understanding of crack opening rates, cement growth rates, and vein-forming mechanisms is highly relevant in the analysis of naturally fractured geothermal reser- voirs (e.g., Dezayes et al. 2010; Haffen et al. 2013) as well as hydrocarbon exploration and production (e.g., Gale et  al. 2004; Hilgers et  al. 2006; Lamarche et  al. 2012; Wuestefeld © The Author(s) 2018. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Prajapati et al. Geotherm Energy (2018) 6:7 Page 2 of 15 et  al. 2016). The presence of open cracks enhances flow rates in geothermal projects and hydrocarbon reservoirs in otherwise impermeable formations. The natural frac - ture growth due to repeated fluid pressurization events (i.e., hydrofracturing) signifies a non-continuous and episodic crack opening that might grow rapidly over human time scales. However, recent petrographic investigations of the fractured sandstone in Cre- taceous Travis Peak formation (Becker et  al. 2010) reveal a prolonged fracture growth driven by gas generation, over a duration of 48 m.y. at the estimated crack opening rates of 16–23 μm/m.y. Due to a large number of physical boundary conditions that are present in different geological settings, there exist several mechanisms that could lead to a specific vein tex - ture. For instance, the early work of Taber (1918) proposed that brou fi s veins are formed due to deposition of mineral in the wall rock-vein interface without any fracturing. This growth mechanism was termed as ‘antitaxial’ by Durney and Ramsay (1973), and was later numerically simulated by Urai et  al. (1991) and Hilgers et  al. (2001), and experi- mentally recreated by Means and Li (2001). On the other hand, Mügge (1928) proposed that, for the fractures in their growth stage, if the cement growth rates are faster than the crack opening rates, the growth competi- tion between neighboring crystals is suppressed and brou fi s textures are formed. In gen - eral, the process of cement growth contemporaneous with fracture opening is known as synkinematic cementation (Laubach 1988). Therefore, any cement that has precipitated at the time of fracture opening, whether a crack-seal texture is present or not, is known as a synkinematic deposit. In nature, the sizes of incremental fracture opening can vary from few microns to few millimeters (Fisher and Brantley 1992). Accounts of accurate incremental opening sizes or so-called ‘gaps’ in the aforementioned range have been fur- ther reported in Laubach et al. (2004b) for quartz and in Hooker et al. (2012) for carbon- ate rocks. During synkinematic cementation, a complicated interplay of the continued crack opening and crystal growth occurs. This interplay determines the morphological attributes of the resulting veins based on the crack opening rates (Mügge 1928; Urai et  al. 1991) and was shown to affect quartz-bearing veins (e.g., Laubach 1988; Lander and Laubach 2015), calcite veins (Hilgers et al. 2001), and carbonate veins in dolostones (Gale et  al. 2010). Carbonate vein textures can also be subject to an interplay of syn- and postkinematic cementation (Ukar and Laubach 2016). The analysis of microtextures using e.g., a cathodoluminescence detector mounted on a scanning electron microscope (SEM-CL) might highlight the different cement generations (e.g., Ukar and Laubach 2016). The natural samples of calcite veins derived from the Jurassic limestones of southern England exhibit four prominent morphologies, as depicted in Fig. 1a–d. Different crys - tallographic orientations of vein crystals are reflected by distinct interference colors. The brou fi s veins (in Fig.  1a) are characterized by parallel grain boundaries occurring as thin fibers of different lengths and high length-to-width ratios. Figure  1b, c depicts the vein microstructures with elongate-blocky crystals manifesting crystallographic preferred ori- entations at the final growth stage, which indicates the presence of growth competition during vein formation. Elongate-blocky veins often exhibit partially open morphologies Prajapati et al. Geotherm Energy (2018) 6:7 Page 3 of 15 Fig. 1 Natural samples of calcite veins in Jurassic limestones from southern England exhibiting a fibrous, b elongate-blocky, c partially open, and d euhedral morphologies. The wall rock in fibrous veins is not contained in the sample, but is located towards the base of the image in a. Different interference colors of vein crystals depict different crystallographic orientations. All images are obtained from transmitted microscopy and crossed polarizers [one of the partially open morphologies is also called rind morphologies by Lander and Laubach (2015)]. Here, euhedral crystal faces border some porosity within the partially open fracture (Fig. 1c). Moreover, there exist veins that exhibit euhedral crystals (charac- terized by straight edges and sharp corners when visualized along thin-section) in con- tact with the pore space (black) in open fracture, as illustrated in Fig. 1d. In this work, we utilize a multiphase-field model to delineate the underlying mecha - nisms behind the formation of the aforementioned calcite-vein textures based on the boundary conditions of synkinematic cementation in limestone. The phase-field method is a diffuse interface approach widely used for modeling microstructure evolution pro - cesses such as phase transitions, grain growth, among others (see review articles Boet- tinger et  al. 2002; Moelans et  al. 2008; Nestler and Choudhury 2011; Johannes et  al. 2016). In contrast with the traditional front-tracking methods (see Urai et al. 1991; Bons 2001; Hilgers et  al. 2001; Nollet et  al. 2005), the phase-field approach does not require explicit tracking of the interfaces as the motion of interfaces is inherently captured in the governing equations. This property makes it a computationally efficient and robust approach for the treatment of moving boundary problems such as cement growth in open fractures or porous rocks, even in three dimensions (e.g., Ankit et al. 2013, 2015a, b; Wendler et al. 2016; Prajapati et al. 2017). The present article is organized as follows: In "Methods " section, we describe the phase-field model, the numerical setup and the algorithm for the geometric shift condi - tion. We then present the simulation results of bitaxial calcite cementation in fractures for different crack opening rates and highlight their implications on the resulting vein textures and the histories of fracture porosity in "Results and discussion" section. Fur- ther, the congruence of our numerically simulated vein textures with the natural samples is discussed. In "Conclusion and outlook" section, we conclude the article by recapitulat- ing the important findings and propose directions for further studies. Prajapati et al. Geotherm Energy (2018) 6:7 Page 4 of 15 Methods Multiphase‑field model A thermodynamically consistent multiphase-field model Nestler et  al. (2005) is employed to address the phenomenon of anisotropic crystal growth in calcite veins. The model has been previously applied to address the growth of minerals such as potash alum (Ankit et  al. 2013), quartz (Ankit et  al. 2013, 2015a, b; Wendler et  al. 2016), and calcite (Prajapati et al. 2017). Therefore, for the sake of completeness, we briefly reiterate the model equations. For a detailed discussion and numerical implementation of equa- tions and process-specific model extensions, the interested readers are referred to the above-mentioned articles. We consider a vector of phase-field parameters φ(x, t) =[φ (x, t), . . . , φ (x, t)] , where 1 N φ ∈[0, 1] ∀ α ∈{1, . . . ,N } describes the presence of each calcite grain α at position x and time t. The grain-grain or grain-liquid interface is characterized by a diffuse region of finite width. In the interface region shared between grain α and liquid/other grains, the phase-field parameter φ varies smoothly from 0 outside, to 1 inside the grain. The location and shape of the grain are defined by the isoline with φ = 0.5 . At each spatial point in the system with computational domain  , a summation constraint φ = 1 α=1 is satisfied. The total free energy of the system is formulated as F(φ,∇φ) = εa(φ,∇φ) + w(φ) + f (φ) dΩ , (1) where ε represents a length scale parameter that controls the interfacial-width. εa(φ,∇φ) , w(φ) , and f (φ) are the gradient energy density, the potential energy density, and the bulk free energy density of the system, respectively. The interface region is described by the gradient energy density, given by N ,N 2 2 εa(φ, ∇φ) = ε γ a (φ, ∇φ)|q | , αβ αβ αβ (2) α,β=1 (α<β) where γ denotes the interfacial energy density between α and β grains, and q αβ αβ (= φ ∇φ − φ ∇φ ) represents the generalized gradient vector pointing in the outward α β β α normal direction of the α–β interface. The choice of function a (φ,∇φ) determines the αβ type of interface energy anisotropy ( a = 1 for the isotropic case). Although, calcite αβ exhibits a wide variety of growth habits including acute to obtuse rhombohedral, sca- lenohedral, prismatic, among several others, we limit the present discussion to a com- mon rhombohedral calcite geometry (see Fig. 2a). In the present study, we assume that the crystal facets are developed due to strong faceted anisotropy in the interface energy. Therefore, a piece-wise function of the form αβ a (φ, ∇φ) = max · η αβ (3) k,αβ 1≤k≤n |q | αβ αβ is utilized for incorporating such faceted anisotropy in the interface energy formula- tion. Here, {η | k = 1, . . . , n } denotes the position vectors corresponding to the n k,αβ αβ αβ Prajapati et al. Geotherm Energy (2018) 6:7 Page 5 of 15 Right−handed Symmetric Left−handed Diffuse 1 interface liquid 90° −90° Top views φ grain grain φ =0.5 grain t=0 t=1505 a Initial spherical grain Simulated shape Rhombohedral calcite geometry Fig. 2 a Symmetric and asymmetric projections of the rhombohedral calcite geometry. b Validation of interface energy anisotropy formulation corresponding to the asymmetric left-handed 2D projection. A spherical grain evolves to its equilibrium shape while retaining a constant volume. The final crystal shares a diffuse interface of finite width with the surrounding liquid. The shape and size of the crystal are defined by the isoline corresponding to φ = 0.5 grain Table 1 Values of phase-field model parameters used in simulations Model parameters Value Grid size x 1 Time-step width t 0.07 Length scale parameter ε 6.5 Interfacial energy density γ , i.e., γ and γ 1 αβ g,g g,l Higher order parameter γ 10 αβδ Kinetic coefficient for grain–grain interface τ 100 g,g Kinetic coefficient for grain–liquid interface τ 1 g,l vertices of the Wulff shape of crystal α dispersed in bulk phase β . The evolution equa - tions of the phase-fields read ∂φ ∂a(φ, ∇φ) ∂a(φ, ∇φ) 1 ∂w(φ) ∂f (φ) τε =ε ∇· − − − − (4) ∂t ∂∇φ ∂φ ε ∂φ ∂φ α α α α for α = 1, . . . ,N . Here, τ is the interfacial kinetic coefficient and  is the Lagrange parameter to ensure the summation constraint φ = 1 . The set of partial differential equations in Eq. α=1 (4) is solved using forward Euler scheme for time derivative and second-order accurate central difference scheme for the spatial derivatives. The numerical solution method is efficiently implemented in a parallel simulation framework for crystal evolution, PACE3D (see Hötzer et al. 2018), which is written in C language. The values of model parameters used in the present work are listed in Table 1. Modeling calcite cements The present work comprises a two-dimensional treatment of the complex process of synkinematic cementation with continuous fracture growth. Therefore, we obtain the following 2D projections: right-handed, symmetric and left-handed (see Fig.  2a), that were also highlighted in our previous article pertaining to three-dimensional mode- ling of cement overgrowth in porous limestones (Prajapati et al. 2017). During bitaxial Prajapati et al. Geotherm Energy (2018) 6:7 Page 6 of 15 cement growth in fractures, the effect of asymmetricity can be significant (Ankit et  al. 2015b). Hence, we chose the left-handed projection of the calcite, in order to account for the influence of asymmetricity. The chosen position vectors for the vertices of crystal shape are as follows: η = [0, ±0.75] η = [1, 0.25]= −η 1,2 3 4 resulting in an aspect ratio of 1.33 of the 2D crystal. Figure  2b depicts the simulated left-handed 2D equilibrium shape of calcite from an initial spherical grain embedded in liquid phase, using the volume preservation technique (Nestler et al. 2008). For sim- ulating bitaxial cement growth in progressive fracture-opening, the numerical setup is described in "Numerical setup" section. Numerical setup The initial grain structure in a 2D computational domain of dimensions 1000x × 1000x is generated using a voronoi diagram based on a random set of points, see Fig. 3a. The gray region below the lower crack surface is rendered stationary and is termed as barrier phase. In the very first step, a fracture of finite width is intro - duced resulting in two crack surfaces with randomly oriented grains as shown in Fig. 3b. The fracture is completely filled with liquid phase (in yellow), which is assumed to be supersaturated with respect to calcite. This assumption is realized by applying a constant driving force ( f ) for the anisotropic crystal growth. In order to repress the influence of the opposing curvature-dependent forces on the crystal growth (see Plapp (2012) for fundamentals of phase-field), a sufficiently high value of f (=− 0.15) is chosen for all the simulations. The defined boundary conditions for the numerical studies are shown in Fig.  3c. The geometric shift boundary condition applied at the lower crack surface is described in "Geometric shift condition" section. The definition of crystal orientations is depicted in Fig. 3d. Geometric shift condition In order to simulate an opening crack, the previous unitaxial crack opening algorithm utilized by Ankit et  al. (2013) was amended for bitaxial growth in syntaxial veins. The present so-called geometric shift algorithm relocates the entire set of grains on the lower surface by a given vertical displacement after a prescribed number of time steps and fills a b c Isolate d 1000 Δ x Liquid Voronoi crack Geometric shift Periodic Periodic 50 Δ x surfaces condition with grains Initial crack introduced orientation Initial Barrier crack width −90 Fig. 3 a Dimensions of the computational domain consisting of a layer of polycrystalline grains (in RGB colors) generated using voronoi diagram and the stationary barrier phase (in gray). b An initial crack filled with liquid (in yellow) is introduced. c The applied boundary condition at the edges of computational domain. d Definition of crystallographic orientation of calcite grains 1000 Δ x Prajapati et al. Geotherm Energy (2018) 6:7 Page 7 of 15 Table 2 Values of  the  prescribed initial crack width and  the  crack opening rates for  three different cases Cases Initial crack width Crack opening rate (COR) Slow opening rate 100x �x/1000�t Moderate opening rate 100x �x/200�t Fast opening rate 100x �x/100�t the new empty cells with liquid phase. This essentially implies that, in the scenario of crystal-bridging or vein closure, the crack always opens along the boundaries between the grains that lie on mutually opposite crack surfaces. Moreover, the phase-field evolu - tion is not computed in the barrier phase, but the barrier is constantly consumed due to the shifting process. Results and discussion Eec ff t of crack opening rate: implications for vein morphologies By employing the generated computational domain with randomly oriented calcite seeds on the crack surfaces (in Fig.  3a), and the geometric shift condition for crack opening (described in "Geometric shift condition" section), simulations were performed for equal initial crack width and different values of crack opening rates, see Table  2 for param- eters. Each simulation was performed until the barrier phase was completely consumed due to the moving lower crack surface along with the attached grains. In particular, three cases of slow, moderate, and fast opening rates were considered as follows: • Slow opening rate: For the applied driving force ( f =− 0.15 ), the prescribed crack opening rate (in Table 2) is smaller than the minimum potential growth. For the cho- sen crack width of 100x (Fig. 4a), the magnitude of initial crack aperture is one order higher than the initial mean grain size of aggregates. This leads to a growth compe - tition in the early stages as illustrated at t = 1190 in Fig.  4b. As soon as the crystals grow and fill the initial crack (at t = 3850 in Fig.  4b), further growth of crystals is almost completely dependent on the opening increments. From this stage onwards, the growth competition arising due to random grain orientations is suppressed. This implies that almost all grains grow equally into the fracture spaces introduced dur- ing small crack opening increments. Different stages of synkinematic cementation are shown in Fig. 4a, b. The resulting vein consists of brou fi s crystals of nearly equal width and parallel to the crack opening displacements (at t = 59605 in Fig. 4b, e). • Moderate opening rate: For the same driving force of grain growth as the previous case, the prescribed crack opening rate (in Table  2) lies between the maximum and minimum potential growth rate. Different stages of synkinematic cementation are shown in Fig. 4c. In the initial stages, due to the availability of sufficient pore space, the growth competition between the neighboring crystals of the same crack surface is dominant. It is noteworthy that since the aspect ratio of the chosen 2D equilibrium shape of calcite is close to unity, the orientation selection is more or less random. Therefore, the fastest growing crystals are decided on the basis of neighboring grain Prajapati et al. Geotherm Energy (2018) 6:7 Page 8 of 15 Early growth competition Suppressed growth due to high initial competition after crack−width first complete sealing Slow crack opening rate Fibers t=1190t=3850 t=59605 e crack opening direction Bridging crystals Moderate crack opening Partially−open rate vein Elongate−blocky t=0.7 t=2660 t=4935 t=11935 f a c Fast crack opening rate Open vein at all stages Euhedral t=1050 t=2450 t=5985 d g Fig. 4 a Computational domain with an initial crack width of 100x . Simulated bitaxial crystal growth for different values of crack opening rates: b Slow crack opening rate �x/1000�t, c moderate crack opening rate �x/200�t, and d fast crack opening rate �x/100�t . The progression is shown at representative simulation time steps. The selected crystals in the resulting vein morphologies are illustrated. e Fibers with almost equal width. f Elongate-blocky crystals with increasing width. g Crystals forming euhedral flat facets due to availability of sufficient pore space for crystal growth interactions. As the cementation progresses, the crystals begin to form bridges (at t = 2660 in Fig.  4c), so that the further growth is constrained by the crystals grow- ing from the opposite crack surface resulting in partially open veins (at t = 4935 in Fig.  4c). At later stages, the pore space is entirely filled with bitaxial crystal growth leading to elongate-blocky vein textures with increasing width along the length of the crystals (at t = 11935 in Fig. 4c, f ). • Fast opening rates: In this case, the crack opening rates are higher than the maxi- mum potential growth rate. Under such conditions, the growth competition between the neighboring crystals is present and the growth is not constrained by the crys- tals growing from the opposite fracture surface. The veins are open at all stages of time and the resulting crystals are euhedral in nature with flat facets and sharp cor - ners, as depicted at t = 5985 in Fig. 4d, g. The texture resembles rind morphologies described in Lander and Laubach (2015). In general, the aforementioned simulated results for calcite cement growth in fractures, depending upon fracture opening rate and crystal orientations are in good agreement with natural and experimental results, presented previously for quartz cements in fractures (Lander and Laubach 2015). Eec ff t of crack opening rate: implications on porosity evolution The phase-field simulations (in "Effect of crack opening rate: implications for vein morphologies" section) provide valuable insights into the processes controlling the Prajapati et al. Geotherm Energy (2018) 6:7 Page 9 of 15 formation of fibrous, partially open, elongate-blocky, and euhedral veins based on the boundary conditions of synkinematic cementation. Apart from recreating the mech- anisms behind various vein morphologies, the numerical simulations can provide information on histories of fracture porosity as a result of synkinematic cementation for different crack opening rates, see Fig.  5. The porosity-time behavior results from the interplay of crack opening rate and the anisotropic growth of crystals. For the present study, constant crack opening rates are applied for all cases, thereby implying a linear porosity-time behavior in the absence of anisotropic crystal growth (plotted in same colored dashed lines as the fracture sealing curves in Fig. 5). The anisotropic growth of crystals adds the non-linearity in the time evolution of porosity for the three crack opening rates (slow, moderate and fast). We analyze each case one by one as follows: • Slow opening rate (crack opening rate = �x/1000�t ): The porosity-time curve monotonically decreases with respect to time. Due to the early growth competition (described in "Effect of crack opening rate: implications for vein morphologies " sec- tion), the initial decline in the porosity is slow but accelerates as the cementation proceeds. After the first complete sealing (somewhere close to t = 3850 as shown in Fig. 4b), the crack reopening occurs along the interfaces between the crystals grow- ing from opposite sides of the vein. The subsequent porosity evolution is completely dependent on the crack opening increment. The temporal evolution of porosity enters a periodic oscillatory regime due to the occurrence of alternate events of crack opening and vein filling, as illustrated in the zoomed inset picture of Fig. 5. • Moderate opening rate (crack opening rate = �x/200�t ): The porosity initially increases with time owing to the early growth competition (that suppresses crystal growth) and a higher opening rate. The growth competition declines due to occlu - sion of the growth of unfavorably oriented crystals, resulting in an increase in the net crystal growth rates. These increased growth rates lead to a maximum (at t = 1225 ) followed by a monotonic decrease in the porosity-time curve. Due to the higher crack opening rate than the former case, the first complete sealing takes place at later stage in the present case. 0.8 Δx/ 1000Δt Δx/ 200Δt 0.7 Δx/ 100Δt Δx/ 50 Δt t=5985 0.6 t=3010 0.5 0.4 t=59605 0.3 t=0.7 0.2 Periodic 0.131 0.1 t=11935 0 2000 4000 6000 8000 10000 12000 59605 1225 time Fig. 5 Temporal evolution of fracture porosity for different crack opening rates represented by different colored solid lines. The same colored dashed lines depict the porosity-time behavior in the absence of cement growth, signifying constant crack opening rates. The inset picture illustrates the initial and final vein morphologies at the end of simulation porosity Prajapati et al. Geotherm Energy (2018) 6:7 Page 10 of 15 • Fast opening rate (crack opening rates = �x/100�t and �x/50�t ): The evolution of porosity monotonically increases with time. For even larger opening rates (i.e., �x/50�t ), the porosity-time curve gets steeper. It is noteworthy that the porosity- time curve has an upper bound that is prescribed by the crack opening rates (dashed lines in Fig. 5). Therefore, although the porosity increases monotonically, the poros - ity-time curve is concave with respect to time. Eec ff t of initial crack aperture From the previous set of simulations (in "Effect of crack opening rate: Implications for vein morphologies" section), it is evident that the growth competition of different grain orientations is suppressed by the slow crack opening rates leading to fibrous morphol - ogies. However, the early growth competition that arises due to the large initial crack aperture (one order higher than the grain size) results in the consumption of less favora- bly oriented grains in the initial stages. In this section, we investigate the influence of initial crack width on the resulting fibrous morphologies. In particular, we consider two cases of small and large initial cracks. The parameters used for the present study are listed in Table 3. – Small initial crack: In this case, the magnitude of initial crack width and grain size are of same order. Under such conditions, the growth competition is suppressed from the very beginning, leading to retardation in the occlusion of unfavorably oriented grains by the favorable ones. Therefore, more initial grains are preserved as shown at t = 8400 in Fig. 6a. As the cementation proceeds, the unfavorably oriented crystals are gradually occluded after the initial stages, resulting in arrested fibers of reduced length at intermediate time step ( t = 31150 in Fig.  6a). The resulting veins (at t = 66220 in Fig. 6a) possess non-uniform fibrous textures consisting of fibers of different lengths. – Large initial crack: The magnitude of crack width is one order higher than that of grain size. Therefore, the intense initial growth competition (at t = 1890 in Fig.  6b) leads to occlusion of almost all the unfavorably oriented grains in the early stages. The fib - ers formed as a result of growth of favorably oriented crystals during the intermediate stages (at t = 24535 in Fig.  6b) exhibit very little growth competition. Therefore, the resulting veins (at t = 59605 in Fig. 6b) have uniform fibrous morphologies composed of almost all fully grown non-occluded fibers. In the following sections, we refer to these fibers as the ‘survivor’ fibers that have continued to grow from the crack surface to the center of the vein, without being occluded by the neighboring crystals. We plot the temporal evolution of porosity along with the initial and final vein mor - phologies for the above cases in Fig.  7a. For a small initial crack width, the opening Table 3 Values of the prescribed initial crack width and the crack opening rate Cases Initial crack width Crack opening rate Small initial crack 5x �x/1000�t Large initial crack 100x �x/1000�t Prajapati et al. Geotherm Energy (2018) 6:7 Page 11 of 15 t=66220 Survivor fibers Small crack t=31150 t=8400 t=0.7 ICW: 5 Δ x Negligible growth More intial grain competition orientations are preserved Arrested fibers VW: 177 Δ x VW: 500 Δ x VW: 1000 Δ x Almost all survivor fibers Large crack b Initial growth Fewer initial grain orientations are competition preserved ICW: 100 Δ x t=0.7 t=1890 t=24535 ICW: Initial crack width Increasing vein width vein−width (VW) t=59605 Fig. 6 Simulated bitaxial crystal growth for a initial crack width = 5x (small crack) and b initial crack width = 100x (large crack). Progression is shown at representative values of the evolving vein-width ( VW ). For small cracks, the resulting veins (at t = 66220 in a) consist of arrested as well as survivor fibers of different lengths, while large cracks lead to uniform fibrous vein textures (at t = 59605 in b) 0.12 70 Small initial crack Small initial crack crack opening rate width (=5 Δx) Δx) width (=5 Δx/ 1000Δt 0.1 Large initial crack Large initial crack width (=100 Δx) width (=100 Δx) 0.08 t=0.7 t=59605 t=66220 0.06 0.04 Periodic 0.02 0 20 200 300 400 500 600 700 800 900 1000 0 2000 4000 60000 66220 time vein width [ x] Δ a t=0.7 b Fig. 7 a Temporal evolution of fracture porosity for different initial crack widths. The inset pictures illustrate the initial and final vein morphologies at the end of simulation. b Plot of the number of survivor grains with the evolving vein-width for different values of initial crack width increments are of same order, forcing the cement growth to be entirely dependent on the crack opening. Therefore, the porosity evolution enters the periodic regime at an ear - lier stage than for the larger crack opening. As highlighted in Fig. 6, the simulated over- growths for the two cases are compared at same vein-width (VW). In order to analyze the influence of initial growth competition, we plot the number of survivor grains with the evolving vein-width, see Fig. 7b. For a vein-width of 200x , the number of survivor grains is lower for the large ini- tial crack due to the suppression of growth of unfavorably oriented grains in the very beginning. The remaining grains grow with little competition; thereby, the decline in the number of survivor grains is lower as compared to the case of small initial crack. For the small initial crack, the unfavorably oriented grains are not occluded initially and are reduced gradually over the entire process of cementation. It is noteworthy that the resulting number of survivor fibers is similar in both cases. Present simulations provide valuable insights into the growth dynamics of arrested and survivor fibers. We infer that natural samples of non-uniform fibrous veins exhibiting arrested and survivor fibers porosity number of survivor grains Prajapati et al. Geotherm Energy (2018) 6:7 Page 12 of 15 result due to slow opening rates and small crack apertures at all stages of synkinematic cementation. On the other hand, uniformly long fibers in the fibrous veins indicate the presence of large crack apertures at the initial stages of synkinematic cementation, lead- ing to the occlusion of unfavorably oriented crystals in the very beginning. Comparison with natural samples Natural samples of calcite veins showcase a wide variety of vein cement morpholo- gies, which can be linked to the simulation results. All shown samples were observed in Jurassic limestones from southern England (Blue Lias at Kilve and Lyme Regis, hand specimens were collected from loose boulders at the beaches). Euhedral vein-form- ing crystals, characterized by straight edges in contact with open pore space (black in Fig.  8a), are in good agreement with simulated examples, highlighting growth compe- tition at the contact with the host rock and increasing crystal width towards the pore (Fig. 8a). Pore morphologies and occasional irregular domain contacts between bitaxially grown calcite vein cements in partially open veins show a good match between natural and simulated examples (Fig. 8b). This includes the possibility of euhedral crystal facets in small open pores. The elongate-blocky morphology of moderate crack opening rates, resulting in a mostly sealed vein is characterized by irregular contacts between domains grown from opposing sides of the crack (Fig. 8c). This also includes growth competition at the host rock interface of the fracture. The case of slow crack opening rate and small initial crack aperture resulting in non-uniform fibrous calcite cement of varying length matches well with habits observed in natural brou fi s veins (Fig. 8d). Conclusion and outlook While the recent numerical works (e.g., Ankit et al. 2015a, b; Wendler et al. 2016) advo- cate the potential of the phase-field approach in addressing quartz mineral growth in sandstones, the present study in addition to our previous work (Prajapati et  al. 2017) serves as a basis for new calcite cementation models for carbonate rocks based on sound thermodynamic principles. The present study complements the works of Bons (2001) and Hilgers et  al. (2001) (who used front-tracking technique) and adds to them. While the vein morphologies arising due to the interplay of different crack opening rates and Fast crack−opening rate Moderate crack−opening rate Slow crack−opening rate a Euhedral veins Partially−openaElongate−blocky d Fibrous b c suppressed grains Euhedral crystal facets in small open spaces suppressed grains due to growth Euhedral crystals Vein fibers of different lengths competition at host rock interface Fig. 8 Comparison of simulated vein morphologies with natural samples. The sketches of natural samples (in gray scale) are drawn from the natural samples of calcite veins in Jurassic limestones from southern England, shown in Fig. 1. a Euhedral, b partially open, c elongate-blocky and d fibrous veins. The wall rock (gray) in the sketch of fibrous vein is not contained in the sample, but it is towards the base of the image. The simulated veins exhibit strong resemblance with the natural samples suggesting a pronounced role of crack opening rates and initial crack width in the formation of different vein textures Prajapati et al. Geotherm Energy (2018) 6:7 Page 13 of 15 cement growth rates are faithfully reproduced by the phase-field simulations, the present study additionally highlights the formation mechanisms of uniform and non-uniform fibrous veins based on the initial crack apertures. The convincing resemblance between the numerically sealed vein morphologies with natural samples further demonstrates the capabilities of the employed multiphase-field model in elucidating the growth dynam - ics of four prominent calcite-vein textures (fibrous, elongate-blocky, partially open, and euhedral). Furthermore, the model in conjunction with innovative pre- and post-pro- cessing tools developed for the analysis and visualization of numerical data allows to speculate the reduction of fracture porosity over time. Although the present study is limited to 2D, the implementation of the model in 3D is fairly straight forward, albeit with increased computational resources. This provides a distinct edge over the conventional front-tracking methods that are restricted to 2D owing to numerical issues that arise while addressing the process in three dimensions. Although bitaxial crack-seal mechanism is by far the most common vein-forming pro- cess, nevertheless, very few numerical studies have hitherto been conducted, due to the aforementioned reasons. The extension of the present work to three dimensions remains an interesting prospect for future studies that enables the analyses of the fracture net- works and flow behavior in the geothermal reservoir rocks. The present work assumes a constant crack opening and models fracturing as a continuous process. However, the implemented geometric shift condition can directly be utilized to simulate episodic reopening of fractures (orthogonal to the crack surface) in response to repeated pres- surization events, i.e., hydrofracturing. The present work paves the way for the ongoing studies that explore the development of syntaxial textures due to complex fracture open- ing histories and paths such as crack opening obliquely to the vein boundaries or along curved trajectories. Moreover, building upon the previous works of (Ankit et al. 2013, b), the model along with appropriate extensions can be utilized to further explore the controls on the formation of antitaxial and ataxial veins. The present work is an initial treatment of a highly coupled thermo-hydro-mechanical-chemical problem of vein for- mation as a result of cementation, wherein faceted crystal growth mechanisms based on the boundary conditions of synkinematic cementation are simulated. The present phase- field model needs to be extended in order to account for the mechanical forces for crack opening and chemical driving forces for crystal growth that facilitates the differential growth rates on facetted and unfacetted surfaces as suggested by Ukar and Laubach (2016). Provided the availability of physical parameters (e.g., crystal orientations, super- saturation of solute, absolute rates of calcite accumulation) through benchmark cement growth experiments, the model can be calibrated to precisely mimic the experimentally synthesized veins and further explore the controls on fabric development in natural veins. Authors’ contributions NP conducted the numerical studies under the supervision of MS in addition toconceptual discussions with BN, BB, and CH. MS and NP implemented the geometric shift algorithm for the present study. BB and CH provided the micropho- tographs of the natural samples. NP drafted the manuscript. All authors editedthe manuscript. All authors read and approved the final manuscript. Author details Institute for Applied Materials (IAM-CMS), Karlsruhe Institute of Technology, Building 30.48, Strasse am Forum 7, 76131 Karlsruhe, Germany. Institute of Digital Materials Science (IDM), Karlsruhe University of Applied Sciences, Molt- kestr. 30, 76133 Karlsruhe, Germany. Institute of Applied Geosciences (AGW-SGT ), Karlsruhe Institute of Technology, Adenauerring 20a, 76131 Karlsruhe, Germany. Prajapati et al. Geotherm Energy (2018) 6:7 Page 14 of 15 Acknowledgements The authors would like to thank the Helmholtz Association for the financial support through “KIT-Integration initiative in geothermal systems” and within the program “RE-renewable energies.” Further, we acknowledge publication support by Helmholtz Centre for Environmental Research—UFZ; Helmholtz Centre Potsdam—GFZ German Research Centre for Geosciences; and Karlsruhe Institute of Technology for the article processing charges. Competing interests The authors declare that they have no competing interests. Availability of data and materials All relevant data and material are presented in the paper. Consent for publication Not applicable. Ethics approval and consent to participate Not applicable. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Received: 15 March 2018 Accepted: 26 May 2018 References Ankit K, Nestler B, Selzer M, Reichardt M. Phase-field study of grain boundary tracking behavior in crack-seal microstruc- tures. Contrib Miner Petrol. 2013;166(6):1709–23. Ankit K, Selzer M, Hilgers C, Nestler B. Phase-field modeling of fracture cementation processes in 3-D. J Petrol Sci Res. 2015;4(2):79–96. Ankit K, Urai JL, Nestler B. Microstructural evolution in bitaxial crack-seal veins: a phase-field study. J Geophys Res. 2015;120(5):3096–118. Bathurst RG. Carbonate sediments and their diagenesis, vol. 12. Amsterdam: Elsevier; 1972. Becker SP, Eichhubl P, Laubach SE, Reed RM, Lander RH, Bodnar RJ. A 48 my history of fracture opening, temperature, and fluid pressure: Cretaceous Travis Peak Formation. East Texas basin. Bulletin. 2010;122(7–8):1081–93. Boettinger WJ, Warren JA, Beckermann C, Karma A. Phase-field simulation of solidification. Ann Rev Mater Res. 2002;32(1):163–94. Bons PD. Development of crystal morphology during unitaxial growth in a progressively widening vein: I. The numerical model. J Struct Geol. 2001;23(6–7):865–72. Bons PD, Elburg MA, Gomez-Rivas E. A review of the formation of tectonic veins and their microstructures. J Struct Geol. 2012;43:33–62. Cox SF, Etheridge MA. Crack-seal fibre growth mechanisms and their significance in the development of oriented layer silicate microstructures. Tectonophysics. 1983;92(1–3):147–70. Dezayes C, Genter A, Valley B. Structure of the low permeable naturally fractured geothermal reservoir at Soultz. Comptes Rendus Geosci. 2010;342(7–8):517–30. Durney D, Ramsay JG. Incremental strains measured by syntectonic crystal growths. In: De Jong KA, Scholten R, editors. Gravity and tectonics. New York: Wiley; 1973. p. 67–96. Fisher DM, Brantley SL. Models of quartz overgrowth and vein formation: deformation and episodic fluid flow in an ancient subduction zone. J Geophys Res. 1992;97(B13):20043–61. Gale JF, Laubach SE, Marrett RA, Olson JE, Holder J, Reed RM. Predicting and characterizing fractures in dolostone reser- voirs: using the link between diagenesis and fracturing. Geol Soc Lond. 2004;235(1):177–92. Gale JF, Lander RH, Reed RM, Laubach SE. Modeling fracture porosity evolution in dolostone. J Struct Geol. 2010;32(9):1201–11. Girard JP, Sanjuan B, Czernichowski-Lauriol I, Fouillac C. Diagenesis of the Oseberg Sandstone Reservoir (North Sea): an example of integration of core, formation fluid and geochemical modelling studies. AAPG Bulletin. 5(CONF-960527). Gutjahr A, Dabringhaus H, Lacmann R. Studies of the growth and dissolution kinetics of the C aCO polymorphs calcite and aragonite II. The influence of divalent cation additives on the growth and dissolution rates. J Cryst Growth. 1996;158(3):310–5. Haffen S, Géraud Y, Diraison M, Dezayes C. Fluid-flow zones in a geothermal sandstone reservoir: localization from ther - mal conductivity and temperature logs, borehole EPS1(Soultz-sous-Forêts, France). Geothermics. 2013;46:32–41. Herman JS, Lorah MM. Calcite precipitation rates in the field: measurement and prediction for a travertine-depositing stream. Geochimica et Cosmochimica Acta. 1988;52(10):2347–55. Hilgers C, Koehn D, Bons PD, Urai JL. Development of crystal morphology during unitaxial growth in a progressively wid- ening vein: II. Numerical simulations of the evolution of antitaxial fibrous veins. J Struct Geol. 2001;23(6–7):873–85. Hilgers C, Sindern S. Textural and isotopic evidence on the fluid source and transport mechanism of antitaxial fibrous microstructures from the Alps and the Appalachians. Geofluids. 2005;5(4):239–50. Hilgers C, Urai JL. Microstructural observations on natural syntectonic fibrous veins: implications for the growth process. Tectonophysics. 2002;352(3–4):257–74. Prajapati et al. Geotherm Energy (2018) 6:7 Page 15 of 15 Hilgers C, Nollet S, Schonherr J, Urai JL. Paleo-overpressure formation and dissipation in reservoir rocks. Oil Gas Eur Mag. 2006;32(2):68–73. Hilgers C, Kirschner DL, Breton JP, Urai JL. Fracture sealing and fluid overpressures in limestones of the Jabal Akhdar dome. Oman mountains. Geofluids. 2006;6(2):168–84. Hooker JN, Gomez LA, Laubach SE, Gale JFW, Marrett R. Eec ff ts of diagenesis (cement precipitation) during fracture open- ing on fracture aperture-size scaling in carbonate rocks. Geol Soc Lond. 2012;370(1):187–206. Hötzer J, Reiter A, Hierl H, Steinmetz P, Selzer M, Nestler B. The parallel multi-physics phase-field framework Pace3D. J Comput Sci. 2018;26:1–12. Johannes H, Michael K, Philipp S, Britta N. Applications of the phase-field method for the solidification of microstructures in multi-component systems. J Indian Inst Sci. 2016;96(3):235–56. Lamarche J, Lavenu AP, Gauthier BD, Guglielmi Y, Jayet O. Relationships between fracture patterns, geodynamics and mechanical stratigraphy in Carbonates (South-East Basin, France). Tectonophysics. 2012;581:231–45. Lander RH, Laubach SE. Insights into rates of fracture growth and sealing from a model for quartz cementation in frac- tured sandstones. GSA Bull. 2015;127(3–4):516–38. Laubach SE. Subsurface fractures and their relationship to stress history in East Texas Basin sandstone. Tectonophysics. 1988;156(1–2):37–49. Laubach SE, Reed RM, Olson JE, Lander RH, Bonnell LM. Coevolution of crack-seal texture and fracture porosity in sedi- mentary rocks: cathodoluminescence observations of regional fractures. J Struct Geol. 2004;26(5):967–82. Laubach SE, Lander RH, Bonnell LM, Olson JE, Reed RM. Opening histories of fractures in sandstone. Geol Soc Lond. 2004;231(1):1–9. McNamara DD, Lister A, Prior DJ. Calcite sealing in a fractured geothermal reservoir: insights from combined EBSD and chemistry mapping. J Volcanol Geotherm Res. 2016;323:38–52. Means WD, Li T. A laboratory simulation of fibrous veins: some first observations. J Struct Geol. 2001;23(6–7):857–63. Moelans N, Blanpain B, Wollants P. An introduction to phase-field modeling of microstructure evolution. Calphad. 2008;32(2):268–94. Mügge O. Ueber die Entstehung faseriger Minerale und ihrer Aggregationsformen. Neues Jahrbuch für Mineralogie, Geologie und Paläontologie A. 1928;58:303–438. Nestler B, Choudhury A. Phase-field modeling of multi-component systems. Curr Opin Solid State Mater Sci. 2011;15(3):93–105. Nestler B, Garcke H, Stinner B. Multicomponent alloy solidification: phase-field modeling and simulations. Phys Rev E. 2005;71(4):041609. Nestler B, Wendler F, Selzer M, Stinner B, Garcke H. Phase-field model for multiphase systems with preserved volume frac- tions. Phys Rev E. 2008;78(1):011604. Nollet S, Urai JL, Bons PD, Hilgers C. Numerical simulations of polycrystal growth in veins. J Struct Geol. 2005;27(2):217–30. Plapp Mathis. Phase-field models. Multiphase microfluidics: the diffuse interface model. Vienna: Springer; 2012. p. 129–75. Plummer LN, Parkhurst DL, Wigley TML. Critical review of the kinetics of calcite dissolution and precipitation. In: Jenne E. Chemical modeling in aqueous systems; ACS Symposium Series. 1979, p. 537-73. Prajapati N, Selzer M, Nestler B. Computational modeling of calcite cementation in saline limestone aquifers: a phase- field study. Geotherm Energy. 2017;5(1):15. Ramsay JG. The crack-seal mechanism of rock deformation. Nature. 1980;284:135–9. Taber S. The origin of veinlets in the Silurian and Devonian strata of central New York. J Geol. 1918;26(1):56–73. Ukar E, Laubach SE. Syn-and postkinematic cement textures in fractured carbonate rocks: insights from advanced cathodoluminescence imaging. Tectonophysics. 2016;690:190–205. Urai JL, Williams PF, Van Roermund HLM. Kinematics of crystal growth in syntectonic fibrous veins. J Struct Geol. 1991;13(7):823–36. Wendler F, Okamoto A, Blum P. Phase-field modeling of epitaxial growth of polycrystalline quartz veins in hydrothermal experiments. Geofluids. 2016;16(2):211–30. Wuestefeld P, De Medeiros M, Koehrer B, Sibbing D, Kobbelt L, Hilgers C. Automated workflow to derive LIDAR fracture statistics for the DFN modelling of a tight gas sandstone reservoir analog. In: 78th EAGE Conference and Exhibition 2016; 2016.

Journal

Geothermal EnergySpringer Journals

Published: Jun 2, 2018

References

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off