Quartz Cementation in Polycrystalline Sandstone: Insights From Phase‐Field Simulations

Present work investigates the dynamics of polycrystalline quartz cement growth in sandstone using a multiphase‐field model. First, the model parameters corresponding to common reservoir temperature and pressure conditions were determined. A parameter related to growth kinetics was ascertained through undisturbed cement growth simulations to aptly capture the known faceting‐dependent growth behavior of quartz. Unrestricted growth simulations for different grain sizes, number of subgrains, and their crystallographic orientations revealed that (I) the model successfully recovers the tendency of quartz cements to grow at a faster overall rate on a coarse grain as compared to finer one and (II) the impact of crystallographic orientations of individual subgrains in polycrystalline grains on cement volume increases with increasing number of subgrains. For applying the model to realistic multigrain systems, we generated digital grain packs through a systematic procedure. These packs fairly represent natural sandstone in terms of grain shapes, sizes, and depositional porosity. The simulated textures in these packs resemble natural samples in terms of crystal morphologies and pore geometries. The cement growth rate decreases with the increase in fraction of polycrystalline grains, indicating a significant role of mutual hindrance among differently oriented overgrowths originating from different subgrains. Further, the permeabilities computed using fluid‐flow simulations through the progressively cemented packs suggest that monocrystalline packs are more permeable than equally porous polycrystalline ones. The computed permeabilities of the simulated microstructures are consistent with the measured permeabilities, advocating the pack generation process and cementation modeling approach.


Introduction
The reduction of porosity in arenitic sandstone is governed by physical and geochemical processes that include mechanical compaction and cementation (e.g., Paxton et al., 2002;Taylor et al., 2004;Worden & Burley, 2003;Waugh, 1970;Walderhaug, 2000). The intergranular pore space is occluded by authigenic cement phases (such as quartz) resulting in the reduction of storage and flow capacity of sandstone reservoirs and thus reservoir quality (Amthor & Okkerman, 1998;Ajdukiewicz & Lander, 2010;Oelkers et al., 1996;Taylor et al., 2015;Worden et al., 2018). Under elevated temperatures above 75 • C, quartz cement overgrowths may precipitate syntaxially on detrital grains from silica supersaturated fluids (Lander et al., 2008;Walderhaug, 1994aWalderhaug, , 1994b. These cements form in optical continuity with the grains and inherit their crystallographic orientations. Natural textures of quartz cement growing anisotropically in porous sandstones composed of monocrystalline quartz grains were successfully modeled by Prajapati et al. (2018b). However, depending on the source area, many quartzarenites also contain polycrystalline detrital quartz grains consisting of several subdomains of quartz crystals (or "subgrains"). Each subgrain produces different overgrowths in accordance with the crystallographic orientation of the polycrystalline detrital grain and attains euhedral form when sufficient pore space is present (Lander et al., 2008;McBride, 1989;Waugh, 1970).
For example, Figure 1 depicts the photomicrographs of an aeolian sandstone (BQ 1 from Busch et al., 2017). The pigmented hematite dust rims, represented by brownish lines in Figure 1a, separate the detrital grains from the optically continuous cement overgrowths. Dust rims thus provide information about the original rounded detrital grain shapes. When visualized under crossed polarizers in Figure 1b, it is found that detrital grains are composed of monocrystalline and polycrystalline grains. The figure indicates that under identical thermal exposure (i.e., time spent at elevated temperatures) and silica supersaturation of the pore fluid, the amount of quartz cement formed on equal-sized monocrystalline and polycrystalline grains possessing same surface area also differs. Quartz cement growth experiments in hydrothermal reactors (Heald & Renton, 1966) have revealed that monocrystalline grains grow at a faster rate than polycrystalline aggregates and the rate reduces with decreasing size of the subgrains in the latter case. This phenomenon is in accordance with the study pertaining to natural samples of Nugget Sandstone (James et al., 1986), which reported that monocrystalline grains are more likely to have well-developed overgrowths as compared to polycrystalline counterparts. Worden and Morad (2000) suggested that the competitive growth between separate and differently orientated quartz subgrains is responsible for the common observation that overgrowths are less extensively developed on polycrystalline grains than on single crystals. This essentially implies that for polycrystalline detrital quartz grains with identical crystal structure, the quartz overgrowth volume differs based on crystallographic orientations of subgrains and reduces with increasing mutual interference between the outgrowths from neighboring subgrains. Therefore, sandstones composed of larger fractions of polycrystalline and chert grains are expected to have less quartz cement.
It is well known that conventional techniques such as well logging provide only a limited understanding of the quartz cementation process, histories, and future pore space evolution. For providing deeper insights into the process of quartz cement growth under different growth conditions, hydrothermal experiments have played a useful role (Ballman & Laudise, 1963;Cecil & Heald, 1971;Heald & Renton, 1966;Lander et al., 2008;Okamoto & Sekine, 2011;Okamoto et al., 2010;Pittman, 1972;Thomas et al., 1949;Van Praagh, 1947). Several experimental studies have confirmed that quartz growth is fastest along the c axis as compared to the a axis (Van Praagh, 1947;Pittman, 1972;Lander et al., 2008). As a result, elongated quartz crystals are formed under undisturbed growth conditions. Ballman and Laudise (1963) determined the following order of growth rates of the common quartz facets: prismatic m{1010} < pyramidal z{1101} < pyramidal r{0111} < basal plane {0001}. Among pyramidal facets, the slower growth rate of z facets has also been confirmed in the laboratory synthesis (Lander et al., 2008) as well as in the microinfrared spectroscopic analysis (Ihinger & Zink, 2000) of quartz crystals. Moreover, the experiments of Lander et al. (2008) revealed that the pre-euhedral growth in quartz is more rapid than the post-euhedral and found that growth rate along the c axis drops by a factor of 20 after the euhedral faceting is complete. Their results indicated that this tendency of quartz cements is responsible for the overgrowths to grow more rapidly on larger grains compared to the smaller ones.
Numerical approaches, in parallel or as viable alternatives to experimental investigations, can serve as powerful tools for evaluating the relationships between petrophysical properties and the evolving microstructures in sandstone-hosted geothermal and hydrocarbon systems that might be difficult to constrain otherwise. Diagenetic modeling tools such as Exemplar (Lander & Walderhaug, 1999) and Touchstone (Makowitz et al., 2006) performed quartz cementation based on input parameters pertaining to rock properties (e.g., available surface area, grain size distribution, and nucleation discontinuities). The corresponding changes in available intergranular volume due to compaction and authigenic cements or matrix are also modeled over time. Recently, these models have been adjusted to also include the impact of monocrystalline versus polycrystalline substrates and a growth rate reduction after reaching euhedral terminations on the volumes of formed syntaxial quartz cements (Lander et al., 2008). Although these models have been proven as sound basis for predrill prediction (Ajdukiewicz & Lander, 2010), they only take into account the globally averaged 10.1029/2019JB019137 data for modeling cementation. Therefore, the morphological details of the evolving microstructures due to anisotropic cement growth are not rigorously accounted. For a specific rock structure, spatiotemporal computational models that are able to simulate the instantaneous changes locally could certainly be advantageous and may provide more trustworthy predictions. Simulation algorithms based on front tracking methods (e.g., Bons, 2001;Hilgers et al., 2001;Zhang & Adams, 2002) provided the capabilities to describe grains with straight faces. They elucidated the complicated grain boundary interferences and mechanisms of formation of natural cement textures in two dimensions (2-D). Later on, Lander et al. (2008) simulated quartz growth in 2-D grain packs using their in-house solver called Prism2D, which was based on the cellular automaton approach. The simulation parameters for growth were calibrated using their own laboratory experiments along with natural sandstone data sets. Their numerical experiments were able to reproduce the common observation of quartz overgrowths to grow at slower rates on polycrystalline grains as compared to monocrystalline ones. These results suggested that slower post-euhedral growth behavior of quartz is majorly responsible for this observation and the minor impact of mutual hindrance among neighboring overgrowths is also present which predates the euhedral effect. Despite the achievements of all the aforementioned works, their scope was restricted to 2-D owing to modeling complexities and computational costs that arise when accounting the process in full 3-D space. Thus, these 2-D models provide limited insights about the petrophysical properties such as permeability which is determined by the nature of the flow passages in 3-D.
In the past decade, the phase-field method has established itself as a viable alternative for simulating mineral growth in 2-D and 3-D (Ankit et al., 2013;Ankit, Urai, et al., 2015;Ankit, Selzer, et al., 2015;Prajapati et al., 2017Prajapati et al., , 2018aPrajapati et al., , 2018bWendler et al., 2011Wendler et al., , 2016. In contrast to the aforementioned front tracking models that were based on the sharp interface approach, the phase-field models describe interfaces as diffuse regions. The increasing utility of these models is due to the fact that they do not require explicit tracking of interfaces, making it a robust and computationally efficient methodology in treating a broad range of moving boundary problems (see review articles Chen, 2002;Moelans et al., 2008;Nestler & Choudhury, 2011;Qin & Bhadeshia, 2010). Prajapati et al. (2018b) analyzed the impact of cementation on the relationships between different rock properties (i.e., porosity, permeability, grain sizes, and pore size statistics) in grain packs composing of monocrystalline spherical grains. However, the model suffered from three main limitations. First, the euhedral shape was generated using interface energy anisotropy alone, and the kinetics was assumed to be isotropic. But it is observed that different facets of quartz grow at different rates (Ballman & Laudise, 1963) and growth along the c axis is more rapid as compared to the a axis (Pittman, 1972;Thomas et al., 1949). Second, the tendency of quartz cements to grow at different pre-and post-euhedral growth rates (Lander et al., 2008) was also not taken into consideration. To account for the abovementioned tendencies, an anisotropy in the growth kinetics is required. Third, a simple interpolation of the inverse mobilities (i.e., ss and sl in Prajapati et al., 2018b) of the solid-solid and solid-liquid interfaces was employed. This interpolation creates an artificial pinning of the triple junctions (shared between the solid-solid-liquid phases) when the inverse mobilities differ by orders of magnitude. To overcome these limitations, in the present work, we extend our previous model following the work of Wendler et al. (2016) and determine the parameter set at the temperature and pressure corresponding to common reservoir conditions. In a systematic manner, we investigate the dynamics of polycrystalline cement growth in digital grain packs composed of realistic grain geometries, to extract quantitative insights of the petrophysical properties.
The present article is organized as follows: Section 2 elaborates the methods employed for the present study. This includes (I) the laboratory analysis of the sandstone sample and (II) the adapted multiphase-field model for cementation. In the model section, we further elucidate the procedure employed to determine the simulation parameter set for the given physical conditions and the workflow utilized for generating the realistic multigrain digital packs. In section 3, we present the simulation results pertaining to undisturbed growth on single crystals as well as constrained growth in multigrain systems. This includes an initial calibration of a growth-related kinetic parameter through undisturbed growth simulations; qualitative comparison of the simulated textures with natural samples; and analysis of the impact of progressive cementation on the dynamic relationships between cement volumes, porosity, and permeability. Finally, we conclude the article in section 4 by recapitulating the main inferences drawn in the numerical investigation and directions for further work.

Methods
Section 2.1 describes the methods employed for analysis of a natural sample. In section 2.2, we delineate the multiphase-field model adapted for simulating quartz cementation. The procedure to determine the simulation parameter set is outlined in section 2.2.1. In section 2.2.2, we present a systematic procedure of generating digital monocrystalline and polycrystalline grain packs composing of grains of realistic shapes that were utilized for numerical investigation of quartz cementation. Finally, in section 2.3, we discuss the utilized fluid-flow model and the applied boundary conditions for computing the permeability of the numerically cemented microstructures.

Laboratory Analysis of Natural Sample
The sample was extracted from Bowscar Quarry (BQ 1 from Busch et al., 2017) and employed for petrographic and petrophysical analysis, as shown in Figure 1. The detrital composition of the sample was determined by point counting (300 steps) on a grid adjusted to the largest grain size and was taken from Busch et al. (2017). Grain sizes were petrographically determined as the long axis of at least 100 grains per sample on a grid adjusted to the maximum observed grain size in the reference area of the sample selected for this study. Sorting was calculated after Trask (1932). Porosity measurements were conducted using a helium pycnometry (micromeritics AccuPyc II 1340). Permeability measurements were conducted at steady state in a Westphal air permeameter at 1.2-MPa confining pressure.

Multiphase-Field Model
Precipitation of silica and subsequent syntaxial overgrowth cementation in reservoir sandstones is energetically governed by the principle of minimization of free energy. Phase-field models can be formulated in a thermodynamically consistent manner to describe crystal growth in such rocks and an entire range of other problems in materials science including solidification and pattern formation in alloys (see Nestler et al., 2005;Stinner et al., 2004). We adapt the multiphase-field model of Wendler et al. (2016) who simulated the epitaxial growth of quartz veins under the experimental conditions of Okamoto et al. (2010) and Okamoto and Sekine (2011). The model equations are recapitulated for two purposes: (I) to enhance the readability and comprehensibility while highlighting the advantages of this model over the previous one (Prajapati et al., 2018b) and (II) to elaborate the methodology of determining the parameter set for temperature and pressure close to reservoir conditions, following the approach proposed by Wendler et al. (2016). The Helmholtz free energy of the system is formulated as where , c, and T are the independent variables. T is the temperature of the system.
represents the set of molar concentrations of K different species.
] denotes a set of N scalar variables known as phase-fields ∈ [0, 1] ∀ ∈ {1, … , N}, each describing the presence of grain (or subgrain in the case of polycrystalline grains) at position x and time t. The interface between each grain/subgrain and liquid/other grains/other subgrains is characterized by a diffuse region of finite width. In this region, varies smoothly from 1 (inside) to 0 (outside) the grain/subgrain . In the right-hand side of equation (1), ( , c, T), a( , ∇ ), and 1 w( ) represent the bulk free energy density, the gradient energy density, and the potential energy density, respectively. Bulk free energy density is given as an interpolation of free energy densities (c, T) ∀ ∈ {1, … , N} of all the bulk phases where h( ) is the interpolation function, continuous in the interval (0, 1), and satisfies h(0) = 0 and h(1) = 1. The difference in the bulk free energy densities of solid and liquid phases (i.e., Δ sl (c, T)= s (c, T)l (c, T)) generates the driving force of crystallization. We assume a constant driving force (i.e., Δ sl = s − l ). This essentially implies that (i) the modeling assumes crystal growth under isothermal conditions and (ii) the solute attachment kinetics at the solid-liquid interface is much slower than diffusion and advection. Therefore, the fluid supersaturation is constant. The diffuse interfaces are energetically accounted by the 10.1029/2019JB019137 gradient energy density contribution given by where denotes the length scale parameter, related to the interface width. The scalar quantity represents the surface energy density of -interface, and q (= − ) is a measure of the phase-field gradient constructed as antisymmetric product of phase-fields and their gradients. For a binary -interface, q points in the direction perpendicular to the interface. The piece-wise function a cap (q ), which is utilized to formulate a faceted-type anisotropy in the surface energy (Nestler et al., 2005), is elaborated in section 2.2.1.4. We chose a multiobstacle type potential energy density formulated as which enforces a sharp minima to the bulk phases through the Gibbs simplex  = { ∈ R N | ∑ N =1 = 1 and ≥ 0}. The second term (∝ ) in equation (4) 1 suppresses the appearance of third phases in the binary interface region (see Hötzer et al., 2016, for a detailed discussion). The evolution equations for N phase-field parameters are derived from the variational principles such that the Helmholtz free energy monotonically decreases with time and is expressed as where  is the variational derivative of the Helmholtz free energy with respect to the vector composing of phase-fields. represents the Lagrange multiplier which enforces the constraint locally. The term ( , ) controls the interface kinetics and is given by after Wendler et al. (2016). 0 is the kinetic mobility of -interface, and ( ) is an interpolation function formulated as if and are solid phases and liquid > 0, The function a kin (q ) introduces anisotropy in the interface mobility accounting for the direction-dependent attachment kinetics and is elaborated in section 2.2.1.4. Several past works (e.g., Ankit et al., 2013;Ankit, Urai, et al., 2015;Ankit, Selzer, et al., 2015;Prajapati et al., 2017Prajapati et al., , 2018aPrajapati et al., , 2018b) employed a simple interpolation of the inverse mobilities (i.e., −1 = . When the mobilities of two interfaces differ in orders of magnitude, this simple interpolation results in an artificial pinning of triple junctions between the solid-solid-liquid interfaces. By using the mobility formulation (equation (6)) and the modified interpolation (equation (7)), the abovementioned numerical artifact is avoided. Substituting equation (1) in equation (5) results in the system of partial differential equations The Lagrange multiplier is computed as = ( ∑ N =1 rhs )∕N. The system of equations (equation (8)) is solved using finite difference schemes (i.e., forward Euler for temporal derivative and second-order accurate central difference scheme for the spatial derivatives) on a regular grid. The model is implemented in the simulation code, Pace3D (2019), using C language. For details about the general structure of solver, its parallelization, and various optimization techniques developed for the computational feasibility of large-scale simulations, interested readers are referred to Hötzer et al. (2018). The procedure for generating the phase-field model parameters at a temperature and pressure corresponding to common reservoir conditions is elaborated in section 2.2.1. The phase-field parameters utilized for the present quartz growth simulations along with their physical values are given in Table 1.

Simulation Parameters for Quartz Growth
We briefly discuss the method adopted to map the phase-field model parameters to their physical values. We consider the temperature of 150 • C as all the required physical parameters are readily available in literature (Okamoto et al., 2010). Generally, the temperature increases at around 30 • C/km with depth (Arndt, 2011). Thus, a temperature of 150 • C in a geological system is expected at the depth of about 4.5 km, when assuming a surface temperature of 15 • C. The hydrostatic pressure according to Bourgoyne et al. (1991) (equation (4.2b) in chapter 4.1) at such depths is around 44 MPa. Therefore, the abovementioned temperature and pressure are considered in the present work. For a detailed discussion regarding the equations, assumptions, and their reasonability, interested readers are referred to Wendler et al. (2016) and the references therein.

Interfacial Energies
We chose the surface energy of quartz-water interface as sl = 0.36 J/m 2 (Parks, 1984). Here s and l denote the solid and liquid phases, respectively. The grain boundary energy ss is determined using the equation for textural equilibrium or Young's law, given by 2 cos For quartz-water system, the dihedral angle has been found to be linearly dependent on temperature (Holness, 1993), with = 60 • at T = 400 • C and = 80 • at T = 600 • C. Therefore, we obtain the ratio ss ∕ sl = 1.91 (or ss = 0.6876 J/m 2 ) at a temperature of 150 • C.

Driving Force of Crystallization
Quartz crystal growth is governed by the precipitation and dissolution of silica (SiO 2 ) described by the chemical reaction where (aq) indicates the dissolved phase. For the crystallization of 1 mole of quartz from a solution, the change in Gibbs free energy of reaction reads (Rimstidt & Barnes, 1980) as a function of gas constant R, solution temperature T (in Kelvins), and fluid supersaturation index S. The driving force of crystallization is given by (Wendler et al., 2016) where  (12) and (13) are based on the following approximations: • The mechanical work (i.e., pΔV = ΔG − Δ ) is negligible as compared to the magnitude of Gibb's free energy, so that equation (11) may directly be utilized for deriving the expression for crystallization driving force, which in the phase-field model, is based on Helmholtz free energy. • The fluid under consideration has low salinity so that supersaturation index can be approximated as This assumption is applicable for dilute solutions when the fluid behavior can be well estimated with that of an ideal solution.
At a temperature T = 150 • C, the equilibrium concentration of silica is C eq H 4 SiO 4 = 60.8 ppm (table A1, Okamoto et al., 2010). For a pressure of p = 44 MPa, the mechanical work during crystallization reads pΔV are the molar volumes of quartz and aqueous orthosilicic acid solution, respectively. At a supersaturation index of S = 1.2 (or C H 4 SiO 4 = 72.96 ppm), a Gibb's free energy change of ΔG = 641.2 J/mol is obtained. In this case, the mechanical work is comparable to the Gibb's free energy change and needs to be included for calculating the driving force. On taking into account the contribution of mechanical work, a driving force of 293.8 J/m 3 is obtained. In the explicit update of the phase-field evolution equation, there exist an upper bound on the simulation time step size (Δt) to ensure the numerical stability of simulations. The maximum allowable time step decreases with increasing kinetic mobility. At the driving force corresponding to S = 1.2, using the maximum allowable time step size (for kinetic anisotropy strength parameter = 115 that recovers the experimentally observed pre-and post-euhedral quartz growth behavior, as discussed in section 3.2.1), the simulation time required to perform the desired amount of quartz cemenation in multigrain systems is unfeasibly large. To tackle this computational constraint, a higher fluid supersaturation can be chosen that would lead to more rapid cementation. At higher supersaturations, the role of model-inherent curvature effects (Plapp, 2012), which account for the inability of crystals above a certain curvature to grow at a given supersaturation, is suppressed. As a result, the microquartz and chert crystals, that would normally cease to grow or would dissolve at supersaturations close to unity, might also grow at higher supersaturations, thereby influencing the resulting morphologies and porosity in natural sandstones. In the knowledge of authors, the extent of this curvature-dependent impact on the polycrystalline morphologies for different fluid supersaturations has not been reported. Therefore, in the absence of reliable data and for the computational feasibility, we assume a supersaturation index of S = 6 (C H 4 SiO 4 = 364.8 ppm) and consider the formation fluid to obey the ideal solution behavior, for our present analysis. For this (p, c, T) growth condition, a crystallization driving force of about 10 5 J/m 3 is generated. We remark that such high values of supersaturation may occur under supercritical conditions in the lower crust and are less likely to be found in the upper crust. Further, the time scale for quartz growth at S = 6 is reduced by a factor of about 10 when compared to the corresponding value at S = 1.2. We further remark that pore waters during quartz cementation may exhibit high salinities, when the low salinity approximation (equation (14)) will break down. Fluid inclusion analyses of the sandstone samples from the same formation report the water salinity values in the range of 13.5 to 23.8 wt% NaCl equivalent (Turner et al., 1995). Therefore, when accounting for the real solution behavior, appropriate correction terms in the calculation of the driving force (equation (12)) will be needed.

Kinetic Mobilities
For quartz precipitation reaction, the interface velocity v sl of solid-liquid interface is given by (Rimstidt & Barnes, 1980) where k − denotes the reaction rate constant for a system consisting of 1 kg of solute and 1 m 2 of quartz surface area and A s ∕M is the area to mass ratio for a given temperature. At higher supersaturations, when the curvature effects in the phase-field model are negligible compared to the crystallization driving force, the kinetic mobility of solid-liquid interface is given by (Wendler et al., 2016) 0 At temperature T = 150 • C, values of all the quantities in equation (16) are provided in Okamoto et al. (2010). We obtain the kinetic coefficient of the solid-liquid interface 0 sl = 4.65 ×10 −17 m 4 /(J s). In order to A representative section was analyzed for extracting the information on grain shapes. Four different grain geometries (G 1 , G 2 , G 3 , and G 4 ), that resemble real grains (when projected in a plane), were chosen. (c) Grain size statistics were obtained by measuring around 500 grains from thin section. (d) Generated pack by simulating sedimentation, that is, free fall of grains in a square domain followed by rearrangement of particles, finally reaching a steady state. (e) Representative cubic volume (edge length = 425 μm) is extracted from the pack and converted to (f) phase-field simulation domain with 165 distinct grain phases, to be utilized for numerical studies.
ensure immobile grain boundaries, we chose a value of 0 ss = 4.65 ×10 −19 m 4 /(J s) for kinetic mobility of the solid-solid interface.

Capillary and Kinetic Anisotropy Functions
The equilibrium shape of a quartz grain in siliceous solution under negligible crystallization driving force is governed by minimization of surface and interface energy. This interface energy anisotropy can be accounted by incorporating a capillary anisotropy function a cap sl (q sl ) of the form wheren = q sl ∕|q sl | is the unit vector perpendicular to the solid-liquid interface. { cap k | k = 1, … , n cap } represents the set of n cap vertex vectors of the resulting capillary anisotropy shape of solid-liquid interface. The function max k returns the largest argument in braces of equation (17). To incorporate a variable attachment kinetics dependent upon crystallographic directions, a kinetic anisotropy function a kin sl (q sl ) of the form (Wendler et al., 2016) Figure 2c. This implies that increasing the value of increases the growth rate along non-euhedral directions, resulting in a higher pre-euhedral growth rate. But, as soon as faceting is complete, the post-euhedral growth is governed by mobilities of the facets, determined by values at cusps in different directions. Value of also controls the maximum allowable time step size Δt for numerical stability of simulations, due to the explicit update of phase-fields. The chosen set of vertex vectors for kinetic anisotropy describes an elongated crystal with an aspect ratio of about 4.6 ( Figure 2d) and imposes a higher post-euhedral growth rate along c-axis direction than a axis. The kinetic anisotropy strength parameter is determined through simulations of undisturbed quartz cement growth to precisely capture the experimentally reported growth behavior (Lander et al., 2008) and is elucidated in section 3.2.1.

Generation of Realistic Digital Grain Packs
Experimental or field samples can be broadly characterized on the basis of (I) size distribution and (II) geometry of the grains. To characterize the grain geometries from a thin section of the same sample (i.e., BQ 1 from Busch et al., 2017, in Figure 3a), we reconstructed the 2-D grain shapes (revealed by reddish dust rims) from a portion of thin section as shown in Figure 3b. Thereafter, from a repository of 3-D digital grains, four Journal of Geophysical Research: Solid Earth 10.1029/2019JB019137 distinct digital grains were chosen, see Figure 3b. These grains sufficiently represent the grains in the thin section when projected in a 2-D plane. The digital grains were generated using a preprocessing algorithm that creates spheres of different sizes in close vicinity based on random function generator and extracts the inner envelop of the overlapping region followed by an edge smoothing process. This process results in randomly shaped 3-D grains. Further, a representative grain size distribution was derived by measuring the maximum distance between two points lying on the grain boundary in the thin section and was approximated with a lognormal distribution function, as depicted in Figure 3c. From the data obtained in the previous steps, simplified sedimentation process was dynamically simulated using Blender (2018) following the belowmentioned four-step preprocessing algorithm: 1. Implementation of the grain size distribution as a mathematical function. 2. Generation of a given number of particles. 3. Assignment of chosen digital grain geometries to particles following the implemented size distribution. 4. Simulation of sedimentation under the action of gravity.
The particle distribution, that composes of digital grains (with geometries arbitrarily chosen from the four digital geometries), follows the grain size statistics (corresponding to the mathematical distribution function) and was allowed to fall under the action of gravity. The grain sizes and geometries play a crucial role, as the simulated deposition involves rearrangement (via tilting and rotation) of grains on collision with each other, until a steady state is achieved. The generated digital unconsolidated pack (in Figure 3d) consists of monocrystalline grains. For the sake of computational feasibility, a representative cubic volume of size 425 μm × 425 μm ×425 μm was extracted ( Figure 3e) and processed into the phase-field data, where each grain was numerically represented by a distinct phase-field variable. The cubic grain pack in Figure 3f composes of 165 grains and exhibits a porosity of about 38%, which is close to the observed porosities of 39-47% for eolian sediments at the time of deposition (Schenk, 1981). In order to investigate the impact of polycrystallinity on the dynamics of cement growth, a polycrystalline and a mixed pack were generated using numerical preprocessing of monocrystalline grain pack and a voronoi structure, see Figure 4. The polycrystalline pack consists of grains with a mean subgrain size of 25.63 μm. The mixed pack is composed of an equal number of monocrystalline and polycrystalline grains and, therefore, represents an intermediate case between the other two packs based on polycrystalline composition. The three digital grain packs were utilized for numerical studies of quartz cement growth.

Computational Fluid Dynamics: Determination of Permeability
Computational fluid dynamics (CFD) analysis was performed to compute the permeabilites by employing the Stokes' model. The digital pore space data at different stages of the numerically cemented digital packs were extracted for the flow simulations. At the solid-liquid interface of the digital pore space, the no slip boundary condition was applied. A pressure drop was prescribed in z direction of the computational domain. In x and directions, the free slip boundary conditions (i.e., zero normal fluid velocities) were applied. The permeabilities were computed using the Darcy's law. The fluid dynamical model equations and the boundary conditions are discussed in Prajapati et al. (2018b). The relevant simulation parameters for fluid flow along with their physical values are given in Table 1. Parameters for flow simulations are utilized solely for evaluating the permeabilities of progressively cemented microstructures and may not necessarily correspond to the fluid-flow conditions in reservoirs.

Results and Discussions
Section 3.1 presents the petrographic and petrophysical information of the sample extracted from the laboratory analysis. The simulation results are presented in sections 3.2 and 3.3. The phase-field parameter set determined from the analysis in section 2.2.1 is listed in Table 1. Some parameters are given without nondimensional values as they are only used in the calculations of crystallization driving force and kinetic mobility. Parameters used in the fluid-flow simulations are also given in Table 1. Section 3.2 includes the results pertaining to undisturbed quartz cement growth on monocrystalline and polycrystalline aggregates. In section 3.3, we present a qualitative comparison of the simulated textures in digital grain packs with natural samples. Finally, the generated numerical data are analyzed to study the influence of polycrystallinity on dynamic relations between cement volumes, porosity, and permeability.

Rock Properties
The presented sample (BQ 1 in Figure 1) is a Permian aeolian sandstone from the Penrith Sandstone Formation, outcropping in the Vale of Eden, Cumbria, N-England. The chosen sample is taken from a dune foreset in Bowscar Quarry north of Penrith. The sample exhibits pinstripe lamination, a sedimentary feature resulting from deposition of coarser grains and finer grains in distinct laminae (e.g., Fryberger & Schenk, 1988), in eolian dune deposits. The average grain size is 0.349 mm, and the sand is moderately sorted. The sample contains 68.3% quartz grains, 2.7% K-feldspar grains, and 0.3% volcanic rock fragments. The quartz grains in the sample compose of 77.07% monocrystalline grains, 5.37% polycrystalline quartz grains with less than five subgrains, and 17.56% with more than five subgrains. The authigenic minerals are 20.7% quartz cement, 0.3% K-feldspar cement, and 2.3% replacements after K-feldspar (illite and kaolinite). Based on the grain composition, the sandstone is composed of 95.8% quartz, 3.8% K-feldspar, and 0.4% volcanic rock fragments. The presented sample can be classified as a quartzarenite after Folk (1980). Due to the large amount of detrital quartz grains, including monocrystalline and polycrystalline grains, this sandstone is a good natural example for the presented modeling approach. The petrophysically determined porosity is at 19.8%, and the Klinkenberg-corrected permeability is at 1,044 mD.

Simulation Results: Undisturbed Cement Growth
Cement overgrowths on quartz grains in reservoir sandstone seldomly result in complete euhedral forms, due to hindrances from neighboring overgrowths in the intergranular pore space. Numerical experiments of undisturbed growth, that physically correspond to cementation on a quartz grain surrounded by an infinitely large pore space filled with siliceous fluid, can be employed to calibrate and validate the model to capture the well-known growth tendencies of quartz cement. For identical phase-field simulation parameters (given in Table 1), unrestricted growth was numerically simulated for different types and sizes of grains, which are listed in columns 1 and 2 of Table 2, respectively. The physical dimensions of the computational domain for each case are depicted in column 3 of Table 2. Each domain was discretized into a regular grid of size Δx = 1 μm. First, we analyze the influence of kinetic anisotropy strength on growth behavior of quartz cements and determine the strength parameter using the experimentally reported growth data. After this calibration, we utilize the full parameter set, which includes the calibrated value of , to investigate the impact of initial grain size, polycrystallinity, number of subgrains, and their crystal orientations on the growth behavior of quartz cements.

Influence of Kinetic Anisotropy Strength
We consider a monocrystalline quartz grain of diameter 90 μm in a simulation domain of size as given in Table 2. The rest of the domain is filled with liquid phase. For the sake of visualization, the quartz grain is depicted in RGB colors, while the liquid phase is represented as transparent, see Figure 5a at t = 0. Using the phase-field simulation parameters given in Table 1 and the set of vertex vectors for anisotropy from Wendler et al. (2016), simulations were performed for different values of the kinetic anisotropy strength parameter . We plot the temporal evolution of velocity along the c axis in Figure 5a. It is observed that velocity along the c axis exhibits maximum value in the initial stages (within 10 days) and its magnitude increases as the value of increases. After a certain growth stage, a constant velocity of 0.00063 mm/day along c axis is attained in all the simulations. The pictures above the plot in Figure 5a depict different stages of quartz cement growth for = 115. As growth progresses, facets begin to appear (at t = 3.45 days), and velocity along the c axis decreases. As soon as faceting is complete (after 20.7 days), velocity along c axis becomes constant. In synthetic quartz precipitation experiments (Lander et al., 2008), it has been found that quartz growth rate along c axis drops by a factor of about 20 after complete faceting occurs. This tendency is recovered for = 115 for which the ratio of initial velocity of the non-euhedral c-face to post-euhedral velocity along c axis is found to be around 20. Predicted growth rates along different facets/directions for = 115 at present growth conditions (p, c, T) are presented in Table 3. It is noteworthy that the predicted growth rate of the reference z facet in simulation matches well with the theoretical value calculated from equation (15). The chosen set of vertex vectors for kinetic anisotropy, that describes a high aspect ratio prismatic quartz crystal, results in highest post-euhedral growth rate along the c axis. This is in qualitative agreement with several experimental studies (e.g., Ballman & Laudise, 1963;Pittman, 1972;Van Praagh, 1947) which confirm that growth of quartz is most rapid along the c axis. As expected, an elongated quartz crystal which is about 0.37 mm long (along c axis) and 0.136 mm wide (along a axis) is obtained after 179.5 days. Further, as the incorporated kinetic anisotropy also accounts for nonequivalent pyramidal r and z facets, at the later stages these facets are clearly visible due to the slower growth rate of z facet than r facet. This tendency of slowly growing z facets has been reported in the early work of Ballman and Laudise (1963) and later confirmed in microinfrared spectroscopic studies of quartz crystals (Ihinger & Zink, 2000) as well as in laboratory synthesis of quartz (Lander et al., 2008). Therefore, for = 115, the simulated growth behavior, both pre-and post-euhedral, is in qualitative agreement with existing literature (Ballman & Laudise, 1963;Ihinger & Zink, 2000;Lander et al., 2008;Pittman, 1972;Van Praagh, 1947).
Unlike free growth, in reservoir sandstones the pore space is limited, and overgrowths might not have sufficient distance to grow before reaching euhedral forms. In this case, an appropriate choice of is crucial as it will influence the pre-euhedral cement volumes. In order to analyze this influence, we plot the temporal evolution of overgrowth cement volume for different values of , as depicted in Figure 5b. The cement growth rate, and thus the overgrowth cement volume, increases with increasing values of as expected. This impact may further increase in multigrain systems where several grains are simultaneously cemented. Hence, we choose the value of = 115 for rest of simulations in this work. The chosen numerical time step Δt in Table 1 is the maximum allowable value for = 115 for ensuring numerically stable simulations. Lander et al. (2008), in their hydrothermal quartz precipitation experiments on different-sized grains, found that quartz cementation occurs at a faster rate on a larger grain compared to a smaller one. The overgrowths growing epitaxially on larger grains require higher amount of cement volumes and must grow farther to attain their euhedral forms. Therefore, the overgrowths of larger grains have the advantage of longer periods of growth on rapidly growing non-euhedral surfaces, leading to this size dependency of quartz cement growth rates. Figure 6a depicts the isosurfaces of the simulated quartz cement growth on a small (grain size = 32 μm) and a large (grain size = 130 μm) monocrystalline spherical grain (domain sizes are given in Table 2). The isosurfaces after an equal amount of time are shown in the same color for both grains. Further, we plot the temporal evolution of absolute and relative overgrowth cement volumes for different grain sizes as shown in Figures 6b and 6c, respectively. The relative overgrowth cement volume is defined as overgrowth cement volume normalized with initial grain volume. We observe the following:

Influence of Initial Grain Size
1. After 3.4 days, the smaller grain already attains flat facets whereas the larger one requires more time before the faceting is complete. 2. As cementation progresses, the gap between overgrowths on smaller and larger grain also increases, as illustrated by colored arrows in Figure 6a. 3. The quartz cement growth rate and cement volume increase with increasing grain size, see Figure 6b. 4. The relative cement growth rate and volume decrease with increasing grain size, as depicted in Figure 6c.
This occurs because with decreasing grain size the surface area per unit volume increases at a rate which overshoots the rate decline due to earlier faceting of smaller grains.
Due to the incorporated kinetic anisotropy formulation with provision for an adjustable pre-euhedral growth rate, the present model is able to recover the size-dependent growth tendencies of quartz cements on non-euhedral grains, consistent with previous experiments (Lander et al., 2008).

Influence of Polycrystallinity of Quartz Grains
In order to investigate the role of c-axis orientations of subgrains on polycrystalline overgrowth cement volumes, two sets of quartz cement growth simulations were performed. The two sets pertain to spherical polycrystalline aggregates composing of two and eight equal subgrains, respectively, which are then compared with simulated overgrowth on a same-sized monocrystalline grain.

Impact of Subgrain Orientation: Monocrystalline Versus Bicrystalline Grains
We consider a monocrystalline grain and the following three cases of c-axis orientations in a bicrystalline grain composed of same-sized subgrains: 1. Orientation difference of 180 • : When the c axes of subgrains are antiparallel to each other. This system physically corresponds to the case of parallel c axis but misaligned a axes resulting in a common solid-solid interface. In this case, the overgrowths are expected to be least hindered. 2. Orientation difference of 135 • : When a certain level of mutual interference is expected to be present. 3. Orientation difference of 90 • : When maximum mutual interference is expected to occur.
Figures 7a-7d depict temporally evolving isosurfaces for growth of monocrystalline and bicrystalline aggregates at representative growth stages. For all cases, the time lapsed at presented growth stages is shown in Figure 7d. Further, we plot the temporal evolution of overgrowth cement volume for all the cases. The cement volume on monocrystalline and bicrystalline grain with antiparallel c axes (Figure 7b) is observed to be nearly equal, as also reflected in the overgrowth cement volume versus time plot in Figure 7e. In the latter case, local curvatures exist near the solid-solid-liquid triple junctions due to differing solid-solid and solid-liquid interfacial energies (i.e., ss ∕ sl = 1.91), illustrated in inset of Figure 7b. However, their impact on overgrowth volumes is negligible. On comparison of cement volumes for the three bicrystalline cases, a small but systematic difference in the overgrowth volumes is observed. In the initial stages, cement volume is lower for the aggregate with an orientation difference of 90 • compared to the one with a difference of 135 • . The pattern, however, reverses at later stages, as illustrated in inset of Figure 7e. This behavior can be rationalized based on two opposing factors: (a) the level of mutual interference which suppresses cement growth and (b) the area of unhindered pyramidal surfaces as they grow faster than prismatic facets. Cement volume plots indicate that the impact of mutual hindrance is predominant in the initial stages, leading to retardation in growth when the c axes of subgrains are perpendicular to each other. Further, for the presented set, it can be inferred that in the third case when c-axis orientation differ by 135 • , in spite of low initial hindrance, at later stages the formed faster growing facets could not attain large enough surface area so as to compete with the free growth of the last case with more number of faster growing facets. The present set of simulations indicate that an interplay of the aforementioned factors governs the quartz cement volume on bicrystalline grains and might even result in larger overgrowth volumes on bicrystalline aggregates than single crystal grains when the impact of unhindered and more rapidly growing surfaces dominates the mutual hindrance.

Impact of Subgrain Orientation: Monocrystalline Versus Polycrystalline Grain
The following three cases of subgrain orientation states of polycrystalline grains were considered: 1. Favorable: When all subgrains can grow freely with least possible hindrance from neighboring subgrains. 2. Unfavorable: When all subgrains experience maximum hindrance from the neighbors. 3. Random: Any configuration lying between favorable and unfavorable orientation states.
Cement volume on a grain with unfavorably oriented subgrains is found to be significantly reduced as compared to favorable configuration, see Figure 8. Essentially, the available surface area for overgrowth formation is highest for the favorably oriented case and lowest for the unfavorably oriented state, which can be directly inferred from the simulated growth ((I) and (III) in Figure 8). In the favorably oriented configuration, minimal hindrances and larger area of faster growing pyramidal surfaces collectively account for highest quartz cement growth rates and vice versa for the unfavorable configuration. Unlike in bicrystalline grains, the difference in cement growth rates and overgrowth volume between favorable and unfavorable configurations is larger. Depending upon mutual hindrance and area of unhindered faster growing pyramidal surfaces, cement volume for all the intermediate orientation states including the presented case of randomly oriented subgrains would lie between the bounds of favorable and unfavorable cases, as illustrated by the grey region in Figure 8. It is noteworthy that, although monocrystalline grains have higher pre-euhedral growth rates than individual overgrowths on relatively small subgrains of the same-sized polycrystalline aggregates, the predicted net cement growth rate of the favorably oriented polycrystalline grain is higher. This is due to a larger number of unhindered faster growing pyramidal facets in the latter that collectively dominate the free growth of monocrystalline grain. The result indicates that, despite of influence of rapid pre-euhedral growth on a larger single crystal than smaller subgrains of a polycrystal, relative orientations of subgrains, which directly relates to the level of mutual hindrances among the overgrowths, might be a major factor influencing the cement volumes. Further, a comparative analysis of simulated undisturbed quartz overgrowth on bicrystalline and polycrystalline grains reveals that the maximum possible difference in overgrowth volumes arising due to different subgrain orientations increases with increasing number of subgrains.

Simulation Results: Cementation in Multigrain Systems
In contrast to undisturbed growth, cementation in natural sandstone occurs in a complicated manner that involves interference of overgrowths originating from randomly oriented adjacent grains and additional interactions of the subgrain overgrowths in polycrystalline sandstones. Under identical growth conditions (numerical parameters listed in Table 1), syntaxial overgrowth cementation was simulated in three digital grain packs, namely, monocrystalline, polycrystalline, and mixed packs, that were generated using the procedure described in section 2.2.2. Figures 9a-9i illustrate the simulated microstructures, visualized along a 2-D plane at three different stages of growth. Different grains are colored according to an orientation colormap, which is based on the c-axis tilt or axial tilt of crystals with respect to the positive z axis, as defined in Figure 9l. In the mixed pack, the crystal orientations of grains were inherited from the corresponding grains of the monocrystalline and polycrystalline packs. We analyze the crystal morphologies and resulting pore geometries of numerically cemented microstructures at an intermediate stage (t = 20.7 days). It is observed that monocrystalline grains develop flat facets surrounding the pores, as shown in Figure 9j  and mixed packs resemble the BQ 1 and BQ 2 samples from Busch et al. (2017), respectively, in terms of crystal morphologies, pore geometries, and remaining pore space, see Figures 9m and 9n. Furthermore, grains growing from adjacent planes of the thin section are also recreated in 2-D plane of simulated overgrowths in the mixed pack, which cannot be accounted in similar 2-D simulations.

Implications on Cement Volumes and Porosity
Utilizing the data sets generated through simulations in section 3.3, we plot the temporal evolution of overgrowth cement volume and porosity for the three packs, see Figures 10a and 10b. The amount of cement volume and the accompanying porosity loss is observed to be inversely dependent on polycrystallinity of the packs. For identical grain structure, polycrystalline sandstones have a higher number of differently oriented crystals as compared to monocrystalline ones, while monocrystalline overgrowths are only hindered by those growing on neighboring grains. Polycrystalline overgrowths, however, have to additionally compete with overgrowths of adjoining subgrains within the same aggregate. Therefore, the growth inhibiting impact of mutual interference is expected to be higher in polycrystalline packs than monocrystalline counterparts. The faster pre-euhedral growth tendency of quartz cements may further account for higher overgrowth volumes on larger monocrystalline grains as compared to individual outgrowths on smaller subgrains that already experience higher mutual hindrances. Even if polycrystalline overgrowths have larger total area of unhindered faster growing pyramidal facets during the post-euhedral growth stages, nevertheless these facets will have limited space to grow in multigrain systems to surpass the impact of growth inhibiting factors. Therefore, monocrystalline and polycrystalline packs exhibit highest and lowest net quartz cement growth and porosity reduction rates, respectively. These results are in qualitative agreement with the experimentally quartz-cemented monocrystalline and polycrystalline sands in hydrothermal experiments (Heald & Renton, 1966). It can be logically deduced that in all the possible intermediate grain packs that can be generated by varying the fraction of polycrystalline grains while keeping a constant particle distribution, including the considered mixed pack, the cement volume and porosity would lie in the regime defined by curves of monocrystalline and polycrystalline packs, as illustrated by grey region in Figures 10a and 10b. It is noteworthy that for all packs, the evolution of cement volumes monotonically increases with time. Moreover, all the cement volume versus time curves are concave. The numerical studies of Lander et al. (2008) reported similar trends in the temporal behavior of quartz cements in monocrystalline packs of Figure 11. Plots of permeability over (a) time and (b) porosity (permeability in semilog scale) for monocrystalline, polycrystalline, and mixed packs. At same porosity, permeability of monocrystalline grain pack is higher compared to polycrystalline counterpart. Permeability of mixed pack lies between corresponding values for the other two packs. (c) The fluid velocities along applied pressure drop through the digital packs at 20-22.5% porosity are clearly seen to be highest in monocrystalline pack and lowest in polycrystalline counterpart, as illustrated by the velocity streamlines. different grain sizes. The present study, in addition, further reveals that sandstones composing of polycrystalline grains are expected to exhibit similar kinetic behavior, albeit with lower cement volumes than their monocrystalline counterparts.

Implications on Porosity-Permeability Relationship
We investigate the impact of cementation on the porosity-permeability relationships in sandstones composed of polycrystalline aggregates. From the generated data sets of the three grain packs in section 3.3, digital pore space was extracted at selected time steps. Computational fluid dynamics analysis of the digital pore space was performed in order to compute permeabilities of the progressively cemented microstructures. From the obtained fluid velocities in the direction of applied pressure drop, permeabilities of progressively cemented microstructures were computed using the Darcy's law. We plot the temporal evolution of permeability for the three packs in Figure 11a. The rate of permeability loss is inversely related to the fraction of polycrystalline grains pack, consistent with the predicted trends for porosity loss in Figure 10b. Figure 11b depicts the plot of permeability as a function of porosity for the three packs, where permeability is plotted in logarithmic scale. Simulated cement growth in section 3.3 reduces the porosity (initial value ≈38%, for all cases) in all digital packs to different final values. The permeabilities vary from 7,523.65 mD for uncemented packs to lowest values of about 10 mD. We analyze the permeabilities of the numerically cemented packs at 20-22.5% porosities which are close to the porosity of the natural sample. The experimentally measured values of porosity and permeabilities of the sample are 19.8% and 1,044 mD, respectively. The sample is composed of 22.93% polycrystalline grains and the rest are monocrystalline. At similar composition, the numerical results indicate that the permeabilities are expected to lie in the range of 921.9 mD (value for mixed pack with 50% each) and 1,221.5 mD (monocrystalline pack) for similar porosity values. Walderhaug et al. (2012) proposed a modified Kozeny equation for the estimation of permeability in sandstones with no intergranular clay based on the data fitting analysis of 66 Norwegian continental shelf samples. Their modified Kozeny equation reads where k denotes the permeability in millidarcy (mD), P is the porosity in percent (%), and D represents the grain size in millimeter (mm). For our present sample with porosity of 19.8% and average grain size of 0.349 mm, equation (19) yields a permeability of 1,176 mD. Although the modified Kozeny equation does not take into account the polycrystallinity of grains, it is noteworthy that our sample permeability is well estimated with this equation. These results show agreement between permeabilities predicted by simulations, experiment, and empirical modified Kozeny relation, advocating the efficacy of the presented method of generation of artificial digital grain packs and the employed multiphase-field model for cement growth. The present study further unveils the impact of polycrystallinity on the flow behavior of sandstone reservoirs. At the same porosity, polycrystalline sandstones exhibit lower permeabilities than their monocrystalline counterparts, as indicated by the porosity-permeability curves in Figure 11b. The same inference can be drawn from the fluid velocity fields, as shown in Figure 11c. For equal volume of pore space, the permeability of a structure is dependent on geometrical characteristics like surface area, tortuosity, and connectivity of the flow passages. In polycrystalline packs, a higher number of randomly oriented overgrowths influence the connected pore geometries. The permeability-porosity trends seen in Figure 11b are a result of difference in the connected pore geometries of cemented monocrystalline and polycrystalline packs.

Concluding Remarks
The present work showcases a systematic and thermodynamically consistent approach of constraining the process of syntaxial overgrowth cementation in reservoir sandstones that exhibit a diverse range of grain shapes, sizes, and polycrystalline structure at microscale. First, the phase-field parameter set for temperature and pressure matching common reservoir conditions was determined. The silica supersaturation of the fluid, however, was set considerably higher to expedite the quartz cement growth and reduce the simulation time, which was otherwise impractically large at supersaturations close to unity. The set of vertex vectors of anisotropy that govern the euhedral facet formation in quartz and their relative growth rates were adopted from Wendler et al. (2016). In the undisturbed growth simulations on monocrystalline grains, the ratio of pre-to post-euhedral growth rate along the c axis of about 20, as reported in hydrothermal experiments (Lander et al., 2008), is recovered for a kinetic anisotropy strength parameter of = 115. The post-euhedral 10.1029/2019JB019137 z-facet velocity measured in the phase-field simulation matches well (∼7% difference) with the theoretical value. The order of post-euhedral growth velocities (lowest to highest) along different facets in simulation is in agreement with the existing literature (Ballman & Laudise, 1963;Lander et al., 2008;Ihinger & Zink, 2000;Pittman, 1972). As the known faster pre-euhedral growth tendency of quartz cements is accounted in the modeling, the simulations are able to successfully capture the influence of grain size. Simulations of unrestricted growth on bicrystalline and polycrystalline grains further indicate that net quartz cement volume is influenced by two opposing factors which are (a) mutual hindrances based on different crystallographic orientations of neighboring grains that suppress the growth and (b) area of unhindered faster growing pyramidal facets which facilitates growth. As the number of subgrains increases, the maximum possible difference in the overgrowth volumes (for varying crystallographic orientations of subgrains) increases.
Simulation results of undisturbed growth on monocrystalline and favorably oriented polycrystalline grains indicate that despite of influence of rapid pre-euhedral growth on a larger single crystal than smaller subgrains of a polycrystalline aggregate, if the relative orientations of subgrains are favorable, the net overgrowth volumes might be higher for polycrystalline grains.
Further, we presented a systematic procedure for the generation of synthetic digital packs that mimic real sand structures in terms of grain shapes, size distribution, and depositional porosity. Simulation results of constrained cement growth in these packs are able to capture several petrophysical and petrographical aspects of reservoir sandstones (composed of monocrystalline and polycrystalline grains) undergoing cementation, which are as follows: 1. The numerically cemented microstructures of monocrystalline and mixed packs resemble natural samples of Busch et al. (2017) in terms of crystal morphologies and pore geometries. 2. As the fraction of polycrystalline grains in sandstones increases, the cement growth rate decreases, coherent with the experimental findings of Heald and Renton (1966). In polycrystalline packs, randomly oriented overgrowths have to compete with outgrowths growing on same as well as neighboring grains.
Our numerical results indicate that growth inhibiting factors which include (a) mutual hindrance and (b) earlier attainment of euhedral forms of smaller subgrains are responsible for slower cementation in polycrystalline sandstone. During post-euhedral growth stages, even if the polycrystalline overgrowths have larger overall area of unhindered faster growing facets, these facets will have a disadvantage of limited intergranular pore space to outstrip the growth inhibiting impact of the aforementioned factors. 3. The temporal evolution of cement volume exhibits a concave nonlinear behavior for all packs. Numerical studies of Lander et al. (2008) predicted concave patterns in the quartz cement versus growth step curves of monocrystalline sandstones. The present work further reveals that similar behavior can be expected for sandstones of varying polycrystalline composition, though at lower cement volumes.
Permeabilities at different stages of the progressively cemented microstructures, as determined from the fluid-flow calculations, reveal that: 1. At same porosity, purely monocrystalline sandstones are more permeable than purely polycrystalline counterparts, which is likely due to complex pore geometries in the latter case. 2. The numerically computed porosity-permeability values of the simulated microstructures are in agreement with the experimentally determined values for the natural sample.
Although the present model considers only quartz cement growth in sandstones composing of quartz grains entirely, we remark that the model is capable of simulating growth of other minerals (e.g., calcite and alum) along with inert grain surfaces (due to, e.g., clay coatings). Given the conditional stability of the explicit time update algorithm (for the phase-fields) and the underlying computational limitations, development of unconditionally stable implicit schemes may relax the constraints on the choice of numerical time step. This would enable to choose lower, and more realistic, fluid supersaturations in the phase-field modeling. Further, to incorporate the impact of fluid salinity in the model, the crystallization driving force needs to be appropriately modified and calibrated (based upon, e.g., experimental or field data), such that the deviation of the fluid from the ideal solution behavior is taken into account. The abovementioned extensions, when effectively implemented in the multiphase-field model, would enable the modeling of quartz cementation in sandstones composed of large population of chert and microquartz crystals at realistic fluid supersaturations. The present work paves the way for more advanced cementation models for the purposes like reservoir quality modeling prediction. In addition to addressing the abovementioned open questions, the future works