Multimineral Diagenetic Forward Modeling for Reservoir Quality Prediction in Complex Siliciclastic Reservoirs

ABSTRACT

Reservoir quality (RQ) prediction models for sandstones in use by the oil and gas industry primarily emulate mechanical compaction, quartz cementation, and cementation by a few select aluminosilicate minerals during burial. The modeled cements are treated on a kinetic basis using independent Arrhenius-style rate equations, and nonmodeled cements are constrained by empirical observations and data from analog rocks. The nonmodeled cements pose a significant modeling challenge in complex lithic and arkosic sandstones that contain carbonate, clay, or zeolite cements when analog data are not available for constraint. Twenty-nine samples covering a broad range of sandstone compositions were selected from four fields in four different sedimentary basins to investigate the potential of using modern reactive transport models (RTM) to simulate a more comprehensive suite of minerals involved in sandstone diagenesis. A commercially available RTM, GWB® X1t, was configured to model chemical diagenesis in a time–temperature burial framework conceptually similar to that employed in standard industry RQ forward models. Strategies were developed to consistently constrain kinetic parameters, fluid compositions, and fluid fluxes with the goal of standardizing these parameters to reduce configuration and calibration time. Results show that RTM using such standardized parameters can reproduce relative timing and volume changes caused by mineral dissolution and precipitation that are accurate within the uncertainty of the petrographic data constraint. The results suggest that incorporation of RTM into current industry models could provide a valuable improvement to siliciclastic RQ prediction in frontier plays with complex mineralogies wherever calibration data are sparse or unavailable.

INTRODUCTION

Numerical reservoir quality (RQ) prediction models for sandstone reservoirs have been in broad use by the oil and gas industry for more than two decades (e.g., Geologica’s Exemplar® by Lander and Walderhaug, 1999; Geocosm’s Touchstone model by Lander et al., 2008; ExxonMobil’s reservoir quality forward model [RQFM] by Awwiller et al.; D. N. Awwiller, 2000, personal communication). These models primarily emulate mechanical compaction and quartz cementation of sandstones during burial (Lander and Walderhaug, 1999), which are the primary mechanisms driving porosity and permeability reduction in many highly quartzose sandstone reservoirs. Early efforts at modeling a more comprehensive set of chemical diagenetic reactions in sandstones were attempted using reactive transport models (RTM) in the late 1980s and early 1990s but were largely unsuccessful for a range of reasons that are outlined in part by Ajdukiewicz and Lander (2010). At about the same time RTM were being tested for applications to sandstone diagenesis, key insights into quartz cementation by Walderhaug (1994a) provided the basis for an effective modeling approach to quartz cementation as a function of time and temperature. The quartz cement models presently in use by industry are based on Walderhaug’s (1994b) deductions and employ an Arrhenius-style kinetic equation with no thermodynamic consideration of quartz saturation beyond the simplifying assumption that dissolved silica is available at all times to form quartz (Walderhaug, 1996; Lander and Walderhaug, 1999; Lander et al., 2008). The parameters used to configure the kinetic quartz models typically are adjusted using petrographic data from offset wells or other appropriate analogs for calibration (Taylor et al., 2010). All of the remaining nonquartz cements are constrained by empirical relations developed from petrographic observations and data. Each nonquartz cement phase is introduced into a simulation based on its relative paragenesis as a qualitative or semiquantitative function of time, temperature, or burial depth. Basic kinetic models for illite precipitation (Lander and Bonnell, 2010), chlorite precipitation from the alteration of volcanic rock fragments (Bonnell et al., 2013), and albite precipitation (Perez and Boles, 2005) are available but require calibration similar to the quartz cement model.

Despite the simplifications to the kinetic treatment of quartz cement, the RQ prediction models presently in use by industry typically provide reliable predictions for porosity whenever reasonable burial histories and appropriate rock calibration data are available. A significant gap in the current modeling approach exists, however, in the absence of offset well data for calibration, and this is especially true in the case of more feldspathic or lithic sandstone reservoirs. Feldspars and many types of rock fragments are chemically unstable and are subject to extensive chemical alteration during burial. Diagenetic clay cements, carbonate cements, grain replacements, and secondary porosity commonly result from such chemical alterations, and they can have equal or greater impact on porosity modification than quartz cement. Thus, for arkosic and lithic sandstones in frontier play areas that lack empirical constraint from offset well data, the standard industry models cannot be reliably configured or calibrated to provide good porosity predictions. To address this gap, nonquartz cement predictions from chemical models that are more firmly grounded in first principles and genetic behavior rather than empirical relations are needed.

A genetic approach must consider the fluid compositions, fluid fluxes, and all of the minerals involved in any particular chemical diagenetic system. The complex interactions between aqueous fluids and rock-forming minerals in sedimentary formations can be modeled using established numerical codes known as RTM. The underlying theory behind RTM is well established (e.g., Raffensperger, 1996; Steefel and MacQuarrie, 1996; Bethke, 1997, 2008; Xu and Pruess, 1998; Steefel et al., 2005) and involves solving the coupled equations of nonisothermal fluid flow, multicomponent solute transport (advection, dispersion, and diffusion), and chemical reactions. Historically, RTM has been applied by industry primarily to address diagenesis and RQ problems in carbonate reservoirs (see summary by Xiao and Jones, 2018). In contrast, there are no comparable industry applications of RTM to general RQ predictions in siliciclastic reservoirs, with the exception of a few studies that address individual diagenetic events tied to hydrothermal or other discrete fluid flow events affecting reservoirs (e.g., Lee and Bethke, 1994; Le Gallo et al., 1998; Bolton et al., 1999; Thyne, 2001; Xiao and Jones, 2006; Xiao et al., 2007; Geloni et al., 2015). The primary reasons for their lack of application to general siliciclastic RQ predictions are that (1) none of the commercially available RTM explicitly couple physical compaction, which is a major control on the RQ of sandstones; (2) challenges exist around the configuration of models to produce mineral assemblages that honor the underlying principles from physical chemistry and actual observations from rock data (i.e., calibration); and (3) there are significant uncertainties in setting the boundary conditions, especially for fluid compositions and fluid fluxes.

Figure 1. Sample subsets selected for chemical diagenetic modeling. Samples in set A are quartz arenites to quartz-rich subarkoses. Set B samples range from arkoses to subarkoses. Set C samples range from sublitharenites to litharenites with the lithic fraction dominated by detrital carbonate grains. Set D ranges from feldspathic litharenites to lithic arkoses with the lithic fraction dominated by basic to felsic volcanic rock fragments (per Folk, 1980). F = feldspar; Q = quartz; R = rock fragments.

The fact that RTM can, in principle, accommodate all of the reactions of interest for nonquartz cements and secondary porosity in complex siliciclastic reservoir rocks suggests that there is an opportunity for improved predictive capabilities in siliciclastic diagenetic systems through integration of RTM with current standard RQ modeling technology, with each filling key gaps found in the competing technology’s capabilities. This paper largely focuses on effective configuration and calibration of RTM for the purpose of RQ modeling in siliciclastic reservoirs, with the goal of exploring how well existing RTM can replicate the scope and timing of diagenetic mineral assemblages across a broad range of sandstone compositions while at the same time honoring what is understood about the underlying physical chemistry and hydrogeology of these complex diagenetic systems. If RTM are able to successfully model a wide range of sandstone chemical diagenesis, then there is potential value to incorporating RTM into industry RQ models for improved predictive modeling in mineralogically complex sandstone plays.

SAMPLES AND MODELS

Samples

Four sets of sandstone samples with widely different compositions were chosen for this investigation (Figure 1; Tables 1, 2<). Set A is a quartz arenite to subarkose deposited in a marginal marine environment of deposition containing abundant detrital glauconite and large volumes of quartz cement. Set B is a subarkose to arkose deposited in a deep-water environment of deposition that has undergone substantial feldspar dissolution and associated production of clay cements and grain replacements. Set C ranges from quartz arenites to litharenites that are dominated by detrital carbonate clasts and were deposited in a fluvial–deltaic environment of deposition. Set C also contains trace to pervasive early and late carbonate cements. Set D ranges from quartz-poor feldspathic litharenties to quartz-rich lithic arkoses with the lithic fraction dominated by volcanic rock fragments, and they were deposited in a marginal marine, fluvial–deltaic environment of deposition. Set D contains abundant quartz, carbonate and clay cements, and grain replacements. Texture in the sample sets ranges from very fine to coarse sand, with most being in the very fine and medium sand size range (Figure 2). Mineralogical compositions in sets B and C vary as a function of grain size.

Geological Model

Most RTM use one-dimensional (1-D), two-dimensional, or three-dimensional static geological models defined by fixed grids that do not change in external physical dimensions through a simulation. Typical sandstone RQ forward models used by industry represent a different physical framework in which a representative volume of sediment deposited at the surface is moved through a time-temperature–effective stress framework in a stepwise fashion to represent a typical burial pathway for a sediment package. Along this burial path, the bulk volume of the representative volume is reduced by compaction at the expense of intergranular porosity, and intergranular pore space may be filled with mineral cements depending on the user-defined paragenesis or kinetic models in use.

Figure 2. Textural characteristics of model sets A, B, C, and D. Initial reactive surface areas in the kinetic models (Table 8) are based on grain size considerations and are subsequently reduced by pore-filling cements and grain coatings.

Modeling compaction plus cementation currently is impractical to do with existing RTM because of their fixed-volume frameworks. Coupling compaction to RTM involves picking up the results from each burial step, remeshing with reduced cell volumes, and running the model using this process recursively to represent a complete burial pathway. These types of simulations attempt to couple thermal, hydraulic, mechanical, and chemical (THMC) processes. The only published case in which a THMC-like model has been applied to siliciclastic diagenesis is by Geloni et al. (2015), in which a simple compaction model was coupled to a modified version of TOUGHREACT (Xu and Pruess, 1998).

In addition to the challenges in coupling mechanical compaction to RTM, another hurdle is implementing multistep burial histories in which each burial segment requires constraint of specific thermal, hydraulic, and chemical boundary conditions encountered along a burial history pathway in a sedimentary basin. The requirement of picking up results from one modeling segment and handing it off to the next offers a similar set of challenges as encountered when attempting to couple mechanical compaction with RTM. Geochemists Workbench or GWB® (Bethke, 2008) has recently added a new capability in version 11.0 within GWB’s X1t (a 1-D RTM) that allows one to impose a more complex temperature, fluid flux, and fluid composition history that can be used to approximate a series of changing boundary conditions during burial. The improvements provide for more seamless transitions between each respective modeling segment. It was recognized that the updated model, ignoring compaction, could be employed in a simple one-cell modeling approach incremented through a thermal–hydraulic–chemical framework (Figure 3) to investigate how well RTM might be able to model chemical diagenesis in sandstones over a compete burial pathway. The results presented here examine strategies for setting boundary conditions and test the ability of X1t to reproduce the observed assemblages, volumes, and timing of diagenetic minerals across a diverse range of starting sandstone lithologies during burial.

Figure 3. One-cell modeling approach using X1t (a one-dimensional reactive transport model). The cell is equivalent to the reference volume used in standard industry reservoir quality models and, in a similar fashion, is incremented through a burial history that is used to constrain the thermal–hydraulic–chemical framework relevant to reactive transport model boundary conditions for each burial segment from burial interval 1 to burial interval n. Compaction is not being accounted for in this modeling exercise.

Reactive Transport Model

The theory and modeling approach behind GWB X1t is long established and well described by Bethke (2008). Basic model configuration requires definition of the initial fluid composition and temperature, inlet fluid compositions and temperatures (Table 3), fluid fluxes through time (Table 4), initial model sandstone mineralogy (Table 5), all minerals that are treated on an equilibrium or kinetic basis (Table 6), and the parameters used to configure the kinetic reactions (Tables 7, 8).

xexwzfuqxsxq


It is postulated here that if the genetics of complex diagenetic systems in siliciclastics are systematic, governed by the principles of physical chemistry, and adequately captured by RTM, standard parameterization of these models should be attainable. Standard parameterization is important because it enables predictions that are consistent with actual observations without the need to completely rework the models every time the sediment composition varies from basin to basin or play to play. This is a key test of the usefulness of RTM as applied to siliciclastic RQ prediction, otherwise every potential use would require complete reconfiguration and calibration with each new play that is considered. This would be both time and cost prohibitive for typical industry RQ analysis and prediction efforts.

Initial Sandstone Mineralogy

Starting sandstone mineral compositions representing sediment composition at the time of deposition were estimated for each sample in the calibration data set on the basis of molar volume (Table 5). The starting mineralogies were restricted to various combinations of quartz, potassium feldspar, albite, oligoclase, andesine, pyroxene (Ca-Al pyroxene, hedenbergite, or both), calcite, dolomite, siderite, biotite, glauconite (starting proxy: 50% high Fe-Mg smectite, 50% chamosite), and Ca-smectite (proxy for detrital matrix). The starting quartz volumes were estimated by summing all monocrystalline, polycrystalline, and other quartz grain volumes. Quartz cement volumes were explicitly excluded from the starting composition. Starting feldspar content was based on reported feldspar grain content. Additional intermediate plagioclase was added to the starting composition by back calculating its volume from clay cement and clay grain replacement volumes in the various data sets, assuming conservation of aluminum. The back calculation relied on the reaction stoichiometry of the feldspar to clay (or secondary porosity) reactions, the molar volumes of the reactants and products, and a microporosity factor. Where volcanic rock fragments were present, additional intermediate plagioclase was included based on a fixed percentage (25%) of the total volcanic rock fragment content to represent feldspar microlites. The remaining 75% of the total volcanic rock fragment content represents the glassy component and was modeled as 60% andesine and 40% pyroxene. Starting calcite volumes were based on reported detrital carbonate grain content and fossils. Syndepositional carbonate cements were also included because the current model does not capture the very earliest stages of syndepositional diagenesis. The only clay components included in the starting mineral composition were Ca-smectite as a proxy for detrital matrix and a proxy mixture of 50% high Fe-Mg smectite and 50% chamosite to represent glauconite grains at deposition. Microporosity was accommodated wherever necessary to adjust petrographic mineral volumes to molar volumes that are used in the RTM.

Thermodynamic Database

The standard thermodynamic database in X1t, thermo.tdat, was employed for all modeling. Five mineral solid solutions were added.

  1. Oligoclase (Na0.8Ca0.2Al1.2Si2.8O8);
  2. Andesine (Na0.6Ca0.4Al1.4Si2.6O8);
  3. 3Fe-chlorite (Fe3Mg2Al2Si3O10(OH)8);
  4. 4Fe-chlorite (Fe4MgAl2Si3O10(OH)8); and
  5. 1Fe-biotite (KFeMg2AlSi3O10(OH)2).

Data for these minerals have been calculated as ideal mineral solutions (see Prieto, 2009) using a linear mixing model based on end-member mineral data from thermo.tdat. For oligoclase and andesine, albite and anorthite end-member data were used. For the Fe-chlorite solid solutions, daphnite and clinochlore end-member data were used. The Fe-bearing biotite solid solution was calculated from phlogopite and annite end members.

Kinetic Parameters

The X1t models can accommodate a wide range of rate equations. The work in this paper employed an Arrhenius form of the rate equation to accommodate changes in reaction rate because of changes in temperature and fluid composition during burial. The rate equation included an additional term for feldspar and clay minerals (Table 7) to accommodate the influence of promoting or inhibiting species on the reaction rate. This is important for feldspars in which lower pH conditions cause significant increases in the rate of reaction (see Palandri and Kharaka, 2004). The rate equation has the general form

bltn18066eqd1

Here, rmineral is the reaction rate (mol s−1), AS is the surface area (cm2), A is the Arrhenius pre-exponential factor (mol cm−2 s−1), e is the base of the natural logarithm, EA is the activation energy (J mol−1), R is the gas constant (8.3143 J K−1 mol−1), T is the absolute temperature (K), aj is the activity of the promoting or inhibiting species, and Pj is that species power in the rate law. Q/K is the saturation index (SI).

All kinetic reactions initially were configured in X1t using the experimentally determined kinetic data compiled by Palandri and Kharaka (2004). Disagreement between kinetic data from laboratory experiments and kinetic data obtained from the analysis of natural systems has been discussed by Ajdukiewicz and Lander (2010). They argue that the use of laboratory-derived kinetics limits the predictive capabilities of first-principle geochemical models. Given this potential limitation when applying RTM to RQ forward modeling, the methodology of Lander and Walderhaug (1999) for calibration of kinetic parameters to natural systems has been adopted to calibrate the kinetic parameters used in the X1t models. To accomplish calibration, the Arrhenius pre-exponential factor was adjusted to values that are one to five orders of magnitude smaller than the laboratory values, and the activation energy was left unchanged for most minerals (Table 7). Failure to reduce the reaction rates typically would cause the models to either overpredict the amount of reaction in comparison to the petrographic calibration data or not converge at all. Calibration was acceptable when the results from all simulations largely fell within the error range of the petrographic calibration data. The reduction in the value of the pre-exponential factor here is less than what typically is used with Touchstone or RQFM. Much of this difference can be attributed to the SI term in equation 1. Observe that as the system approaches equilibrium, the value of 1 − SI approaches zero, thus reducing the overall rate of reaction.

Initial surface areas for kinetic reactions in X1t are treated as specific surface areas (Table 8). Quartz and feldspar surface areas were constrained by employing a simple sphere model using the calculation

bltn18066eqd2

where Amineral is the specific surface area (cm2/gmineral) and D is average grain diameter (mm). The starting surface areas for quartz and feldspar are further adjusted (Aadjusted) for early diagenetic cement content and grain coatings using the calculation

bltn18066eqd3

where E is the early carbonate cement fraction and GCF is the estimated grain coating fraction of surface area for framework grains. The surface area factor that reduces effective surface area based on early carbonate cement content, (0.26-E)/0.26, is a simple approximation using a fixed compacted pore volume fraction, 0.26, for a rigid-grain, well-sorted sandstone (Paxton et al., 2002) because of the fact that the models did not accommodate compaction. Early cements in the modeled rocks did not exceed a volume fraction of 0.26. Sandstones that do would require a modification to this calculation. The surface area factor for early cement will require a modified approach in future models that are coupled to physical compaction models. Note also that grain coatings can be mineral specific as is the case in sample set B in which they selectively form on quartz grain surfaces. Reaction further modifies surface area over the progress of a simulation based on a linear relation with mineral volume in the system, which is the default treatment in X1t. Feldspar surface areas are adjusted upward by approximately one order of magnitude to accommodate their propensity for grain fracturing during compaction to capture the impact of that process on reactive surface area (see Ramseyer et al., 1992) where grain coatings or cements may otherwise limit access of fluids to reactive surfaces. All kinetic clay and zeolite minerals were modeled with a 1000 cm2/g value for specific surface area. Carbonates are handled on an equilibrium basis and do not require surface area constraint.

Please log in to read the full article

ABSTRACT

Reservoir quality (RQ) prediction models for sandstones in use by the oil and gas industry primarily emulate mechanical compaction, quartz cementation, and cementation by a few select aluminosilicate minerals during burial. The modeled cements are treated on a kinetic basis using independent Arrhenius-style rate equations, and nonmodeled cements are constrained by empirical observations and data from analog rocks. The nonmodeled cements pose a significant modeling challenge in complex lithic and arkosic sandstones that contain carbonate, clay, or zeolite cements when analog data are not available for constraint. Twenty-nine samples covering a broad range of sandstone compositions were selected from four fields in four different sedimentary basins to investigate the potential of using modern reactive transport models (RTM) to simulate a more comprehensive suite of minerals involved in sandstone diagenesis. A commercially available RTM, GWB® X1t, was configured to model chemical diagenesis in a time–temperature burial framework conceptually similar to that employed in standard industry RQ forward models. Strategies were developed to consistently constrain kinetic parameters, fluid compositions, and fluid fluxes with the goal of standardizing these parameters to reduce configuration and calibration time. Results show that RTM using such standardized parameters can reproduce relative timing and volume changes caused by mineral dissolution and precipitation that are accurate within the uncertainty of the petrographic data constraint. The results suggest that incorporation of RTM into current industry models could provide a valuable improvement to siliciclastic RQ prediction in frontier plays with complex mineralogies wherever calibration data are sparse or unavailable.

INTRODUCTION

Numerical reservoir quality (RQ) prediction models for sandstone reservoirs have been in broad use by the oil and gas industry for more than two decades (e.g., Geologica’s Exemplar® by Lander and Walderhaug, 1999; Geocosm’s Touchstone model by Lander et al., 2008; ExxonMobil’s reservoir quality forward model [RQFM] by Awwiller et al.; D. N. Awwiller, 2000, personal communication). These models primarily emulate mechanical compaction and quartz cementation of sandstones during burial (Lander and Walderhaug, 1999), which are the primary mechanisms driving porosity and permeability reduction in many highly quartzose sandstone reservoirs. Early efforts at modeling a more comprehensive set of chemical diagenetic reactions in sandstones were attempted using reactive transport models (RTM) in the late 1980s and early 1990s but were largely unsuccessful for a range of reasons that are outlined in part by Ajdukiewicz and Lander (2010). At about the same time RTM were being tested for applications to sandstone diagenesis, key insights into quartz cementation by Walderhaug (1994a) provided the basis for an effective modeling approach to quartz cementation as a function of time and temperature. The quartz cement models presently in use by industry are based on Walderhaug’s (1994b) deductions and employ an Arrhenius-style kinetic equation with no thermodynamic consideration of quartz saturation beyond the simplifying assumption that dissolved silica is available at all times to form quartz (Walderhaug, 1996; Lander and Walderhaug, 1999; Lander et al., 2008). The parameters used to configure the kinetic quartz models typically are adjusted using petrographic data from offset wells or other appropriate analogs for calibration (Taylor et al., 2010). All of the remaining nonquartz cements are constrained by empirical relations developed from petrographic observations and data. Each nonquartz cement phase is introduced into a simulation based on its relative paragenesis as a qualitative or semiquantitative function of time, temperature, or burial depth. Basic kinetic models for illite precipitation (Lander and Bonnell, 2010), chlorite precipitation from the alteration of volcanic rock fragments (Bonnell et al., 2013), and albite precipitation (Perez and Boles, 2005) are available but require calibration similar to the quartz cement model.

Despite the simplifications to the kinetic treatment of quartz cement, the RQ prediction models presently in use by industry typically provide reliable predictions for porosity whenever reasonable burial histories and appropriate rock calibration data are available. A significant gap in the current modeling approach exists, however, in the absence of offset well data for calibration, and this is especially true in the case of more feldspathic or lithic sandstone reservoirs. Feldspars and many types of rock fragments are chemically unstable and are subject to extensive chemical alteration during burial. Diagenetic clay cements, carbonate cements, grain replacements, and secondary porosity commonly result from such chemical alterations, and they can have equal or greater impact on porosity modification than quartz cement. Thus, for arkosic and lithic sandstones in frontier play areas that lack empirical constraint from offset well data, the standard industry models cannot be reliably configured or calibrated to provide good porosity predictions. To address this gap, nonquartz cement predictions from chemical models that are more firmly grounded in first principles and genetic behavior rather than empirical relations are needed.

A genetic approach must consider the fluid compositions, fluid fluxes, and all of the minerals involved in any particular chemical diagenetic system. The complex interactions between aqueous fluids and rock-forming minerals in sedimentary formations can be modeled using established numerical codes known as RTM. The underlying theory behind RTM is well established (e.g., Raffensperger, 1996; Steefel and MacQuarrie, 1996; Bethke, 1997, 2008; Xu and Pruess, 1998; Steefel et al., 2005) and involves solving the coupled equations of nonisothermal fluid flow, multicomponent solute transport (advection, dispersion, and diffusion), and chemical reactions. Historically, RTM has been applied by industry primarily to address diagenesis and RQ problems in carbonate reservoirs (see summary by Xiao and Jones, 2018). In contrast, there are no comparable industry applications of RTM to general RQ predictions in siliciclastic reservoirs, with the exception of a few studies that address individual diagenetic events tied to hydrothermal or other discrete fluid flow events affecting reservoirs (e.g., Lee and Bethke, 1994; Le Gallo et al., 1998; Bolton et al., 1999; Thyne, 2001; Xiao and Jones, 2006; Xiao et al., 2007; Geloni et al., 2015). The primary reasons for their lack of application to general siliciclastic RQ predictions are that (1) none of the commercially available RTM explicitly couple physical compaction, which is a major control on the RQ of sandstones; (2) challenges exist around the configuration of models to produce mineral assemblages that honor the underlying principles from physical chemistry and actual observations from rock data (i.e., calibration); and (3) there are significant uncertainties in setting the boundary conditions, especially for fluid compositions and fluid fluxes.

Figure 1. Sample subsets selected for chemical diagenetic modeling. Samples in set A are quartz arenites to quartz-rich subarkoses. Set B samples range from arkoses to subarkoses. Set C samples range from sublitharenites to litharenites with the lithic fraction dominated by detrital carbonate grains. Set D ranges from feldspathic litharenites to lithic arkoses with the lithic fraction dominated by basic to felsic volcanic rock fragments (per Folk, 1980). F = feldspar; Q = quartz; R = rock fragments.

The fact that RTM can, in principle, accommodate all of the reactions of interest for nonquartz cements and secondary porosity in complex siliciclastic reservoir rocks suggests that there is an opportunity for improved predictive capabilities in siliciclastic diagenetic systems through integration of RTM with current standard RQ modeling technology, with each filling key gaps found in the competing technology’s capabilities. This paper largely focuses on effective configuration and calibration of RTM for the purpose of RQ modeling in siliciclastic reservoirs, with the goal of exploring how well existing RTM can replicate the scope and timing of diagenetic mineral assemblages across a broad range of sandstone compositions while at the same time honoring what is understood about the underlying physical chemistry and hydrogeology of these complex diagenetic systems. If RTM are able to successfully model a wide range of sandstone chemical diagenesis, then there is potential value to incorporating RTM into industry RQ models for improved predictive modeling in mineralogically complex sandstone plays.

SAMPLES AND MODELS

Samples

Four sets of sandstone samples with widely different compositions were chosen for this investigation (Figure 1; Tables 1, 2<). Set A is a quartz arenite to subarkose deposited in a marginal marine environment of deposition containing abundant detrital glauconite and large volumes of quartz cement. Set B is a subarkose to arkose deposited in a deep-water environment of deposition that has undergone substantial feldspar dissolution and associated production of clay cements and grain replacements. Set C ranges from quartz arenites to litharenites that are dominated by detrital carbonate clasts and were deposited in a fluvial–deltaic environment of deposition. Set C also contains trace to pervasive early and late carbonate cements. Set D ranges from quartz-poor feldspathic litharenties to quartz-rich lithic arkoses with the lithic fraction dominated by volcanic rock fragments, and they were deposited in a marginal marine, fluvial–deltaic environment of deposition. Set D contains abundant quartz, carbonate and clay cements, and grain replacements. Texture in the sample sets ranges from very fine to coarse sand, with most being in the very fine and medium sand size range (Figure 2). Mineralogical compositions in sets B and C vary as a function of grain size.

Geological Model

Most RTM use one-dimensional (1-D), two-dimensional, or three-dimensional static geological models defined by fixed grids that do not change in external physical dimensions through a simulation. Typical sandstone RQ forward models used by industry represent a different physical framework in which a representative volume of sediment deposited at the surface is moved through a time-temperature–effective stress framework in a stepwise fashion to represent a typical burial pathway for a sediment package. Along this burial path, the bulk volume of the representative volume is reduced by compaction at the expense of intergranular porosity, and intergranular pore space may be filled with mineral cements depending on the user-defined paragenesis or kinetic models in use.

Figure 2. Textural characteristics of model sets A, B, C, and D. Initial reactive surface areas in the kinetic models (Table 8) are based on grain size considerations and are subsequently reduced by pore-filling cements and grain coatings.

Modeling compaction plus cementation currently is impractical to do with existing RTM because of their fixed-volume frameworks. Coupling compaction to RTM involves picking up the results from each burial step, remeshing with reduced cell volumes, and running the model using this process recursively to represent a complete burial pathway. These types of simulations attempt to couple thermal, hydraulic, mechanical, and chemical (THMC) processes. The only published case in which a THMC-like model has been applied to siliciclastic diagenesis is by Geloni et al. (2015), in which a simple compaction model was coupled to a modified version of TOUGHREACT (Xu and Pruess, 1998).

In addition to the challenges in coupling mechanical compaction to RTM, another hurdle is implementing multistep burial histories in which each burial segment requires constraint of specific thermal, hydraulic, and chemical boundary conditions encountered along a burial history pathway in a sedimentary basin. The requirement of picking up results from one modeling segment and handing it off to the next offers a similar set of challenges as encountered when attempting to couple mechanical compaction with RTM. Geochemists Workbench or GWB® (Bethke, 2008) has recently added a new capability in version 11.0 within GWB’s X1t (a 1-D RTM) that allows one to impose a more complex temperature, fluid flux, and fluid composition history that can be used to approximate a series of changing boundary conditions during burial. The improvements provide for more seamless transitions between each respective modeling segment. It was recognized that the updated model, ignoring compaction, could be employed in a simple one-cell modeling approach incremented through a thermal–hydraulic–chemical framework (Figure 3) to investigate how well RTM might be able to model chemical diagenesis in sandstones over a compete burial pathway. The results presented here examine strategies for setting boundary conditions and test the ability of X1t to reproduce the observed assemblages, volumes, and timing of diagenetic minerals across a diverse range of starting sandstone lithologies during burial.

Figure 3. One-cell modeling approach using X1t (a one-dimensional reactive transport model). The cell is equivalent to the reference volume used in standard industry reservoir quality models and, in a similar fashion, is incremented through a burial history that is used to constrain the thermal–hydraulic–chemical framework relevant to reactive transport model boundary conditions for each burial segment from burial interval 1 to burial interval n. Compaction is not being accounted for in this modeling exercise.

Reactive Transport Model

The theory and modeling approach behind GWB X1t is long established and well described by Bethke (2008). Basic model configuration requires definition of the initial fluid composition and temperature, inlet fluid compositions and temperatures (Table 3), fluid fluxes through time (Table 4), initial model sandstone mineralogy (Table 5), all minerals that are treated on an equilibrium or kinetic basis (Table 6), and the parameters used to configure the kinetic reactions (Tables 7, 8).


It is postulated here that if the genetics of complex diagenetic systems in siliciclastics are systematic, governed by the principles of physical chemistry, and adequately captured by RTM, standard parameterization of these models should be attainable. Standard parameterization is important because it enables predictions that are consistent with actual observations without the need to completely rework the models every time the sediment composition varies from basin to basin or play to play. This is a key test of the usefulness of RTM as applied to siliciclastic RQ prediction, otherwise every potential use would require complete reconfiguration and calibration with each new play that is considered. This would be both time and cost prohibitive for typical industry RQ analysis and prediction efforts.

Initial Sandstone Mineralogy

Starting sandstone mineral compositions representing sediment composition at the time of deposition were estimated for each sample in the calibration data set on the basis of molar volume (Table 5). The starting mineralogies were restricted to various combinations of quartz, potassium feldspar, albite, oligoclase, andesine, pyroxene (Ca-Al pyroxene, hedenbergite, or both), calcite, dolomite, siderite, biotite, glauconite (starting proxy: 50% high Fe-Mg smectite, 50% chamosite), and Ca-smectite (proxy for detrital matrix). The starting quartz volumes were estimated by summing all monocrystalline, polycrystalline, and other quartz grain volumes. Quartz cement volumes were explicitly excluded from the starting composition. Starting feldspar content was based on reported feldspar grain content. Additional intermediate plagioclase was added to the starting composition by back calculating its volume from clay cement and clay grain replacement volumes in the various data sets, assuming conservation of aluminum. The back calculation relied on the reaction stoichiometry of the feldspar to clay (or secondary porosity) reactions, the molar volumes of the reactants and products, and a microporosity factor. Where volcanic rock fragments were present, additional intermediate plagioclase was included based on a fixed percentage (25%) of the total volcanic rock fragment content to represent feldspar microlites. The remaining 75% of the total volcanic rock fragment content represents the glassy component and was modeled as 60% andesine and 40% pyroxene. Starting calcite volumes were based on reported detrital carbonate grain content and fossils. Syndepositional carbonate cements were also included because the current model does not capture the very earliest stages of syndepositional diagenesis. The only clay components included in the starting mineral composition were Ca-smectite as a proxy for detrital matrix and a proxy mixture of 50% high Fe-Mg smectite and 50% chamosite to represent glauconite grains at deposition. Microporosity was accommodated wherever necessary to adjust petrographic mineral volumes to molar volumes that are used in the RTM.

Thermodynamic Database

The standard thermodynamic database in X1t, thermo.tdat, was employed for all modeling. Five mineral solid solutions were added.

  1. Oligoclase (Na0.8Ca0.2Al1.2Si2.8O8);
  2. Andesine (Na0.6Ca0.4Al1.4Si2.6O8);
  3. 3Fe-chlorite (Fe3Mg2Al2Si3O10(OH)8);
  4. 4Fe-chlorite (Fe4MgAl2Si3O10(OH)8); and
  5. 1Fe-biotite (KFeMg2AlSi3O10(OH)2).

Data for these minerals have been calculated as ideal mineral solutions (see Prieto, 2009) using a linear mixing model based on end-member mineral data from thermo.tdat. For oligoclase and andesine, albite and anorthite end-member data were used. For the Fe-chlorite solid solutions, daphnite and clinochlore end-member data were used. The Fe-bearing biotite solid solution was calculated from phlogopite and annite end members.

Kinetic Parameters

The X1t models can accommodate a wide range of rate equations. The work in this paper employed an Arrhenius form of the rate equation to accommodate changes in reaction rate because of changes in temperature and fluid composition during burial. The rate equation included an additional term for feldspar and clay minerals (Table 7) to accommodate the influence of promoting or inhibiting species on the reaction rate. This is important for feldspars in which lower pH conditions cause significant increases in the rate of reaction (see Palandri and Kharaka, 2004). The rate equation has the general form

bltn18066eqd1

Here, rmineral is the reaction rate (mol s−1), AS is the surface area (cm2), A is the Arrhenius pre-exponential factor (mol cm−2 s−1), e is the base of the natural logarithm, EA is the activation energy (J mol−1), R is the gas constant (8.3143 J K−1 mol−1), T is the absolute temperature (K), aj is the activity of the promoting or inhibiting species, and Pj is that species power in the rate law. Q/K is the saturation index (SI).

All kinetic reactions initially were configured in X1t using the experimentally determined kinetic data compiled by Palandri and Kharaka (2004). Disagreement between kinetic data from laboratory experiments and kinetic data obtained from the analysis of natural systems has been discussed by Ajdukiewicz and Lander (2010). They argue that the use of laboratory-derived kinetics limits the predictive capabilities of first-principle geochemical models. Given this potential limitation when applying RTM to RQ forward modeling, the methodology of Lander and Walderhaug (1999) for calibration of kinetic parameters to natural systems has been adopted to calibrate the kinetic parameters used in the X1t models. To accomplish calibration, the Arrhenius pre-exponential factor was adjusted to values that are one to five orders of magnitude smaller than the laboratory values, and the activation energy was left unchanged for most minerals (Table 7). Failure to reduce the reaction rates typically would cause the models to either overpredict the amount of reaction in comparison to the petrographic calibration data or not converge at all. Calibration was acceptable when the results from all simulations largely fell within the error range of the petrographic calibration data. The reduction in the value of the pre-exponential factor here is less than what typically is used with Touchstone or RQFM. Much of this difference can be attributed to the SI term in equation 1. Observe that as the system approaches equilibrium, the value of 1 − SI approaches zero, thus reducing the overall rate of reaction.

Initial surface areas for kinetic reactions in X1t are treated as specific surface areas (Table 8). Quartz and feldspar surface areas were constrained by employing a simple sphere model using the calculation

bltn18066eqd2

where Amineral is the specific surface area (cm2/gmineral) and D is average grain diameter (mm). The starting surface areas for quartz and feldspar are further adjusted (Aadjusted) for early diagenetic cement content and grain coatings using the calculation

bltn18066eqd3

where E is the early carbonate cement fraction and GCF is the estimated grain coating fraction of surface area for framework grains. The surface area factor that reduces effective surface area based on early carbonate cement content, (0.26-E)/0.26, is a simple approximation using a fixed compacted pore volume fraction, 0.26, for a rigid-grain, well-sorted sandstone (Paxton et al., 2002) because of the fact that the models did not accommodate compaction. Early cements in the modeled rocks did not exceed a volume fraction of 0.26. Sandstones that do would require a modification to this calculation. The surface area factor for early cement will require a modified approach in future models that are coupled to physical compaction models. Note also that grain coatings can be mineral specific as is the case in sample set B in which they selectively form on quartz grain surfaces. Reaction further modifies surface area over the progress of a simulation based on a linear relation with mineral volume in the system, which is the default treatment in X1t. Feldspar surface areas are adjusted upward by approximately one order of magnitude to accommodate their propensity for grain fracturing during compaction to capture the impact of that process on reactive surface area (see Ramseyer et al., 1992) where grain coatings or cements may otherwise limit access of fluids to reactive surfaces. All kinetic clay and zeolite minerals were modeled with a 1000 cm2/g value for specific surface area. Carbonates are handled on an equilibrium basis and do not require surface area constraint.

Initial plans were to treat all minerals kinetically; however, it was found that convergence became problematic, and the speed of the calculations greatly slowed when trying to treat more than 10 minerals at a time on a kinetic basis. This was addressed mostly by using the approach of Bethke (2008) to select minerals that react rapidly relative to the time span of a calculation for treatment on an equilibrium basis. The use of mineral solid solutions in lieu of modeling two end-member minerals also reduced the number of kinetic reactions. In a few other cases, scenarios were run using kinetic versus equilibrium models to select minerals that could be treated on an equilibrium basis in which it appeared to have minimal impact on the final result.

All models were run with the same calibrated kinetics shown in Table 7 for sample sets A, B, C, and D. The sole exception is for kaolinite in which the kaolinite pre-exponential factor is further reduced by two orders of magnitude when initial unstable iron-bearing minerals in the sandstone exceed approximately 2 vol. %. The data in Table 7 should provide a reasonable starting point when modeling sandstones of similar compositional ranges and reduce or eliminate the time spent calibrating these basic parameters. Surface area and associated rate-controlling factors such as grain coatings and pore-filling cements will still require constraint on a sample-to-sample basis as they do in conventional sandstone RQ forward models.

Conversion of Reactive Transport Model Molar Volumes to Petrographic Volumes

The RTM simulates chemical reactions and the subsequent change in mineral volumes on the basis of molar volume. The calibration data are from optical petrographic analyses and must be corrected for microporosity in the sandstone constituent materials. Microporosity, as it relates to sandstone petrography, includes intergranular and intragranular pores less than approximately 30 μm in diameter that are not readily visible in standard thin sections and thus are not readily captured in standard petrographic analyses. The presence of microporosity manifests in comparisons of petrographically determined visible porosity to helium porosity from routine core analysis in which helium porosity generally exceeds visible petrographic porosity because of helium’s ability to invade micropores. Micropores are readily imaged in the scanning electron microscope (SEM). The microporosity factors used for converting between petrographic and molar data for sample sets A, B, and C were approximated based on default factors used in ExxonMobil’s RQFM. Sample set D was configured using measured microporosity values from a calibrated SEM backscatter contrast technique (R. Klimentidis, 2001, personal communication).

The RTM also do not distinguish pore-filling cements from grain replacements or secondary porosity. Rules based on observations from each of the petrographic data sets plus the corrections for microporosity can be made to proportionate molar volume changes from RTM into petrographic volumes for cements, grain replacements, and secondary porosity; however, only the total diagenetic volumes corrected for microporosity will be considered in this paper. Such rules are not constant and require modification on a project-by-project basis because of variations in the microporosity of various microporous components in sandstones and varying practices by individual sandstone petrographers in classifying grains, cements, replacements, matrix, and porosity types.

Intergranular Volume

Intergranular volume is determined from petrography and is the sum of thin-section intergranular porosity, cement volume, and detrital matrix volume (Paxton et al., 2002). Because compaction cannot presently be accommodated in the models, the present-day intergranular porosity plus cement and matrix volume was used for each sample as a starting point for model porosity. The main impact from a modeling perspective is that intergranular porosity will be significantly underestimated across the early stages of burial and up to the point of full compaction (∼2–3 km [∼6600–9800 ft] of burial; Paxton et al., 2001).

Fluid Flux and Solute Transport

Fluid flux controls advective solute transport. Fluid flux driving mechanisms during burial in sedimentary basins, as discussed by Harrison and Summa (1991), include topographic head-driven groundwater flow during shallow burial, compaction-driven fluid flow during burial to the point of full compaction, and fluid fluxes driven by fluid generation from the illitization of smectitic clays, the maturation of organic matter, or both in later stages of burial. Local fluid fluxes associated with hydrothermal events driven by other mechanisms, such as thermohaline convection, structurally controlled discharge of overpressured fluids, or emplacement of thrust sheets, are also recognized.

All models discussed in this work accommodate head-driven (burial <1 km [<3300 ft]), compaction-driven (burial ∼1–3 km [∼3300–9800 ft]), and fluid generation–driven flow (burial >∼3 km [>∼9800 ft]) regimes at appropriate burial stages based on their respective burial histories. The head-driven and compaction-driven values for fluid flow are constrained using the maximum horizontal fluid flow in these regimes reported by Harrison and Summa (1991) of 58 cm yr−1 (23 in. yr−1) and 6.5 cm yr−1 (2.6 in. yr−1), respectively. Fluid flow from compaction is scaled based on burial rate in which high burial rates plus sedimentary rock not yet at full compaction drive maximum compaction-driven flow. Fluid fluxes driven by fluid generation from the illitization of smectitic clays, the maturation of organic matter, or both were constrained through calibration. In practice, the illitization-related flux was active between 80°C and 130°C in the models. Smaller illitization-related flux values were used in sample set D, which are volcanic-lithic grain rich, in which abundant diagenesis strongly reduces permeability early in a typical burial history. If an illitization-related flux is not accommodated, it is difficult to produce the amounts of quartz cement observed in natural sandstones. Sample set A includes a hydrothermal event driven by compressional tectonics and associated thrusting (Figure 4). Sample set B is from a deep-water fan deposit, and the fluid fluxes were reduced by approximately one order of magnitude in recognition of its more hydraulically closed nature. Sample set C includes a hydrothermal event where deep formation fluids that were in contact with evaporites are expelled up faults into the reservoir during a regional extensional tectonic event (Figure 4).

Figure 4. Temperature versus time (A) and burial depth versus time (B) for the fields from which sample sets A, B, C, and D were obtained. Key hydrothermal events are marked with stars on (A) for field A and field C.

The X1t models accommodate fluid flow either by directly setting a discharge value or by setting a pressure drop from which a flux is calculated based on Darcy’s law and the current model permeability in the time step. The pressure drop approach is the best way to capture reductions in fluid flow as a result of permeability loss from compaction or cementation. For this investigation, however, in which coupling a compaction model was not possible, discharge was set at fixed values for each burial increment. Ranges and rules for the discharge values used in the models are presented in Table 4 and are largely constrained by the fluid velocity ranges presented in Harrison and Summa (1991).

Diffusive transport of solutes also is accommodated in the X1t models. These models all use the 1.0 × 10−6 cm2 s−1 default value.

Fluid Compositions

The RTM simulations are sensitive to fluid composition, especially under high fluid flux conditions. This sensitivity demands that reasonable fluid compositions be chosen for the purposes of modeling. Hanor (1994, 1996) demonstrated that dissolved Ca2+, K+, Mg2+, and Na+ systematically vary as a function of chloride content in a broad range of formation waters. In numerous investigations across many different sedimentary basins, investigators have postulated that metastable thermodynamic buffering of the pore fluids by minerals drive these systematic variations (e.g., Merino, 1975; Nesbitt, 1985; Hanor, 1988, 1994, 1996; Land et al., 1988; Michard and Bastide, 1988; Smith and Ehrenberg, 1989; Hutcheon et al., 1993 and references therein). These trends are used to constrain the major cation composition of the model fluids as a function of salinity for Ca2+, K+, Mg2+, and Na+. Values for dissolved SiO2 and HCO3- content are constrained from temperature-dependent trends derived from the same internal water database. The Al3+ and in situ pH data generally are not available. The method of Palandri and Reed (2001) to correct formation water analyses has been employed to constrain reasonable downhole values for Al3+ and pH, with the assumption that metastable thermodynamic buffering by minerals in the ambient sediment package are governing the values of these parameters. Typically, dissolved aluminum is constrained by mineral phases such as Na-clinoptilolite, smectite clays, illite, kaolinite, chlorite, albite, and K-feldspar, and pH is constrained largely by carbonates and pyrite plus the aluminosilicates that are used to constrain dissolved aluminum. The data set developed here primarily is aimed at sandstone reservoirs that contained dominantly marine pore fluids at the time of deposition, although a mixed fluvial–marine water was accommodated for model set C to model a phase of early diagenesis in a fluvial–marine mixing zone. A broader discussion of the methodology for establishing pore-water compositions will be addressed in a future paper.

Burial History

Burial history information was available from basin modeling work for sample sets A, B, C, and D. The burial histories are summarized in Figure 4. Key hydrothermal events are also annotated on these plots.

Workflow

The basic workflow is presented in Figure 5. Two activities are required to set up X1t with the goal of modeling chemical burial diagenesis. The first is establishing the starting model rock compositions and textures, and the second is establishing the hydrogeological framework in the context of the burial history. The second activity is used to constrain temperature, fluid fluxes, and fluid compositions through time. The model can be run once the boundary conditions are fully defined. X1t output is then converted back into equivalent petrographic volumes. The final output can then be applied to industry RQ forward models to constrain paragenesis and cement volumes, if desired.

Figure 5. The basic workflow to apply X1t (a one-dimensional reactive transport model) to chemical diagenetic predictions in sandstones includes two steps prior to actually configuring and running the model. These two preparatory steps constrain the starting mineral composition and texture for the model rock and constrain the hydrogeological framework. Key considerations in the hydrogeological framework include time–temperature history, controls on surface and burial water compositions, and fluid fluxes. Once this information is in place, the model can be configured and run, and output data can be applied as needed. EOD = environment of deposition; Ss = sandstone.

RESULTS

The goal of the diagenetic modeling effort is to predict a full range of minerals that are known to be relevant controls on RQ. The basic characterization data from thin-section petrography are ground truth for these models, and the success criteria are twofold.

  1. Obtain model realizations for the correct mineral species and volumes across a broad range of sandstone lithologies.
  2. Demonstrate that the results are consistent with the interpreted paragenetic sequences obtained from thin-section petrography.

Predicted Mineral Volumes

The key diagenetic minerals across all four sample sets include calcite, dolomite, siderite, quartz, chlorite, illite–smectite, and kaolinite. Not all mineral phases occur across all four samples sets. Some minerals occur both as cements and as grain replacements. Additionally, secondary porosity from feldspar dissolution is present in all four sample sets, and glauconite dissolution is present in sample set A. No secondary porosity was reported by the petrographer in sample set A, although it is apparent in thin-section images. Figures 6 and 7 are plots of predicted versus observed mineral volumes along with the probable error envelope at the 95.4% confidence level for point count data sets obtained from 250 points (Neilson and Brockman, 1977). Predicted values within this envelope are considered to be accurate within the precision of the point count data that are being used to calibrate the models. The predicted values represent the results of estimated starting compositions (Table 5) that were subsequently modeled by RTM configured for each field and sample by using a standard set of minerals, kinetic parameters, and respective time–temperature–fluid histories. The modeled molar volumes from RTM subsequently were recalculated into petrographic volumes using the microporosity factors and rules that were previously described. Results are presented in Figure 6 for diagenetic calcite, dolomite, siderite, and quartz and in Figure 7 for diagenetic chlorite, kaolinite, illite–smectite, and glauconite. Most of the model results are very good, with the bulk of the data falling within the petrographic error envelope. Note that for the carbonates in sample sets A, B, and D, variations in diagenetic carbonate largely are the result of redistribution of detrital and syndepositional carbonate plus the destruction of the calcic component in plagioclase feldspar. Sample set C contains abundant detrital calcite that is altered to dolomite during early diagenesis and then back to calcite by hydrothermal fluids that also provide additional externally sourced carbonate cement components. The clay predictions are noisier than the other mineral phases. This likely is because of variations in microporosity that are not adequately accounted for by using an average microporosity value. There may also be some challenges in the partitioning of aluminum between the various clay phases.

Figure 6. Predicted versus observed volumes for diagenetic calcite, dolomite, siderite, and quartz appear to be very good with nearly all predictions falling within the petrographic error envelope. Starting syndepositional cement volumes, such as for early calcite or siderite cement, are included in the initial model configurations.

Predicted Paragenesis versus Observed

The paragenetic sequences reported from petrography for significant diagenetic mineral precipitation or dissolution in fields A, B, C, and D, from earliest to latest, are as follows.

A. Siderite and calcite precipitation (syndepositional), feldspar dissolution, chlorite precipitation, illitization of glauconite, quartz precipitation, and minor calcite precipitation.

B. Plagioclase dissolution, illite–smectite precipitation, late kaolinite precipitation, and minor quartz precipitation.

C. Dolomite precipitation as cement and replacements of detrital CaCO3 rock fragments (burial <100 m [<330 ft]), plagioclase dissolution and kaolinite precipitation, calcite precipitation and dedolomitization of early dolomite, and quartz precipitation.

D. Volcanic rock fragment alteration with coeval clinoptilolite, illite–smectite and chlorite grain coat formation and progressive replacement of volcanic rock fragments, calcite precipitation, kaolinite precipitation, albitization of Ca-plagioclase, and quartz and albite precipitation.

Figure 7. Predicted versus observed volumes for diagenetic chlorite, kaolinite, illite–smectite, and glauconite. Glauconite is included with the diagenetic clays because its volume is recalculated from the diagenetic clay data. Results for most of the samples are very good and fall within the petrographic error envelope but have more scatter than the carbonates and quartz. The increased scatter in the predicted diagenetic clay mineral amounts likely results from additional uncertainty in diagenetic volumes introduced by natural variations in microporosity.

Not all phases listed in the paragenetic sequences are present in each individual sample.

The modeled sequence of diagenetic minerals in set A, shown in Figure 8 for sample A-1, is siderite, illite, feldspar dissolution (volume <∼1%, not shown on plot), chlorite, quartz, and calcite (volume <0.4%, not shown on plot). This sequence generally is consistent with the petrographic constraint. The early illite is associated with illitization of the glauconite proxy. The timing in the model is too slow for syndepositional illitization, as discussed by Hower (1961), and too rapid for normal illitization. The kinetics may be off, or the proxy starting minerals for glauconite may not be adequate. The model also captures some subtleties in the paragenesis. Note that iron chlorite replaces siderite late in the paragenesis at circa 15 Ma. This was not mentioned by the original petrographer, but upon re-examination of the thin-section images, it is clear that early sphaerosiderite has been replaced by Fe-chlorite, as shown in Figure 9.

Figure 8. Sample A-1 model paragenesis. The sequence of events from the model are consistent with the sequence reported by the original petrographers working field A (see text). The model also captures some subtleties in the paragenesis. For example, note that iron chlorite replaces siderite late in the paragenesis at circa 15 Ma. This was not noted by the original petrographer but is evident in thin section (see Figure 9). I–S = illite–smectite.

For field B, the modeled sequence of diagenetic alteration is plagioclase dissolution, kaolinite precipitation along with minor illite–smectite, and minor quartz. An example of this sequence is shown in Figure 10 for sample B-3. Rapid burial at 5 Ma drives a late fluid event that results in additional late feldspar dissolution and production of additional quartz and kaolinite cements. This sequence is consistent with the petrography from field B.

Figure 9. Sphaerosiderite preserved in early calcite cement is shown in (A). Fe-chlorite after sphaerosiderite is interpreted in (B) based on the similarity of the replacement’s size, spherical morphology, and distribution on grain surfaces compared with samples containing preserved sphaerosiderite. Although this is a minor event from a reservoir quality perspective, it is important to note that the simulated reaction sequence captured in reactive transport models can make accurate predictions of even subtle aspects of diagenesis when these models are adequately configured.

The sequence of modeled events in field C is dolomite cementation and dolomitization of detrital calcite (close to time of deposition), followed by calcite cement and calcite replacement of dolomite during a hydrothermal event at circa 134 Ma, minor feldspar dissolution plus kaolinite precipitation, and significant quartz cementation. An example of this sequence is provided in Figure 11 for sample C-2 and is consistent with the petrographic observations of samples from field C. Kaolinite and feldspar are not shown because of their small volumes (<0.5%).

Figure 10. The paragenetic sequence for field B, sample B-3 largely conforms to the sequence of events reported in field B petrography (see text). I–S = illite–smectite.

The modeled sequence in field D is dissolution of volcanic rock fragment proxy constituents (hedenbergite, Ca-Al pyroxene, followed by andesine then oligoclase), Fe-chlorite precipitation, Na-clinoptilolite and illite–smectite precipitation, and quartz and albite precipitation. An example of this sequence is shown in Figure 12 for sample D-4 and is consistent with the petrographic observations of samples from field D.

Figure 11. Modeled paragenesis for field C, sample C-2. The major diagenetic events of early dolomite cementation and grain replacement followed by a hydrothermal dedolomitization event are captured as well as the production of minor quartz cement and kaolinite. These results conform well with the sequence of events reported in field C petrography (see text).

The four paragenetic sequences presented above are on single sample models. They are largely consistent with the general paragenetic sequences observed in other samples from each of the fields.

xexwzfuqxsxq Figure 12. Modeled paragenesis for field D, sample D-4. This is a more complex paragenesis than in the previous examples but is consistent with the sequence of events reported in field D petrography (see text). I–S = illite–smectite; VRF = volcanic rock fragment.

DISCUSSION

The results for model sets A, B, C, and D demonstrate that good predictions of mineral dissolution and precipitation, including volumes and timing, can be made for sandstones with widely different starting compositions and complex mineralogies. Accurate predictions of both major and more subtle changes in mineralogy and the timing of their occurrences during diagenesis suggest that RTM effectively captures many of the genetic processes that control chemical diagenesis in sandstones when boundary conditions are adequately constrained. Furthermore, this result was accomplished using an RTM in which all mineral phases are treated on a common basis in terms of their thermodynamic and kinetic characteristics. The value of attaining a standard configuration for the models is the reduced need for complete reconfiguration for every new project. Additionally, the X1t model configuration files for each of these sandstones may be used as preconfigured starting points for sandstones of similar compositions in similar geological settings. A standard set of modeling fluids for marginal to deep-marine sandstones also has been defined and used along with an approach to constraining the fluid flux history of the sedimentary rocks. General rules have also been developed for converting molar volumes used in RTM to equivalent petrographic volumes and for distributing the petrographic volumes into their respective cements and replacements, although details of these rules are specific to each sample set.

Modeling Considerations and Limitations

A common criticism in the past of RTM applied to siliciclastic diagenesis is that they produce minerals not observed in actual siliciclastic rocks. An important step to obtaining realistic model answers that agree with real-world observations is suppressing minerals that are not known to form under typical conditions associated with burial diagenesis. It was found that key minerals to suppress include high-temperature phyllosilicates (excluding biotite and annite) and ordered dolomite. Thermodynamic data for disordered dolomite were used in lieu of plain dolomite or ordered dolomite options in the thermo.tdat data set. During the initial calibration work, it was found that the reported Mg2+ concentration trends in formation waters from siliciclastic reservoirs could not be honored unless this step was taken. In the model results, there is a suggestion that at temperatures somewhere above 120°C, ordered dolomite may be the preferable phase to use, but this has not been investigated at this point. The 25 minerals used in the final simulations are shown in Table 5, and all others have been suppressed.

The lack of compaction feedback on porosity and permeability also is a key concern. The practice of setting a fixed flux rate for individual burial segments produced acceptable answers for the testing that was performed here, but for effective routine use in predictive RQ modeling for sandstones, a coupled compaction model is needed to attain expected dynamic reductions in fluid flux that result from both compaction and cementation processes or, conversely, increases in fluid flux if sufficient secondary porosity is developed.

The need to reduce the kaolinite pre-exponential factor by two orders of magnitude in Fe-bearing sandstones suggests that the current model configuration is not fully capturing the genetics of Fe in the modeled system during simulated diagenesis. Although the current practice of changing the kaolinite rate works, it would be better if the model operated in a fashion that captures the changes in the resultant diagenesis without intervention. This area will need further attention.

Translation between petrographic volumes and molar volumes requires a thorough understanding of the calibration data sets and the way in which grains, cements, replacements, depositional matrix, and porosity types were characterized. In the case of modeling a new compositional suite of sandstones not previously considered, a calibration exercise from analog rocks likely will be required that will also require working through the exercise of recalculating appropriate starting compositions. In the case that the forward modeling exercise to constrain RQ is from a lithological composition and hydrogeological framework that has already been considered and calibrated, the starting sandstone composition can simply be estimated, or more complex tools, such as SandGEM (Heins and Kairo, 2007), could be employed to constrain starting mineralogy and texture. Note that unstable heavy minerals must be considered and included as part of this activity. Translation of RTM output from molar volumes to petrographic volumes will be required in all cases, except when the RTM output is used only to define the timing and suite of diagenetic minerals for configuration of the diagenesis model in conventional industry RQ forward models. Translation is aided in both directions by the use of basic spreadsheets.

Lastly, the hydrogeological system considered here is restricted to marginal- and deep-marine settings with pore fluids mostly derived and evolved from normal seawater. Systems dominated by nonmarine fluids or high-salinity brines will require additional effort to expand the existing fluid library. It is anticipated that the approach used for the marginal marine system can be effectively used to define new model fluids for other environments. The strategy moving forward is to broaden the initial rock and fluid libraries established in the current work so that the results from previous efforts can be used to rapidly set up new models for routine project support.

CONCLUSIONS

The results from this research have shown that it is possible to use calibrated RTM to perform full mineralogical predictions of chemical diagenesis across a reasonably broad range of sandstone compositions in marine and marginal marine settings. The results of the modeling have been tested and calibrated using four reservoir sandstones with good results. These models, even in their current state of development, can be used to constrain the paragenesis models in standard industry RQ forward models. Ultimately, RTM integrated with a comprehensive compaction model could provide a potential step change in RQ prediction technology for frontier plays with complex mineralogies in which current technology is less useful in the absence of calibration and more fully genetic models are desirable.

REFERENCES CITED

Ajdukiewicz, J. M., and R. H. Lander, 2010, Sandstone reservoir quality prediction: The state of the art: AAPG Bulletin, v. 94, no. 8, p. 1083–1091, doi:10.1306/intro060110.

Bethke, C. M., 1997, Modelling transport in reacting geochemical systems: Comptes Rendus de l’Académie des Sciences, v. 324, p. 513–528.

Bethke, C. M., 2008, Geochemical and biogeochemical reaction modeling, 2nd ed.: Cambridge, United Kingdom, Cambridge University Press, 564 p.

Bolton, E. W., A. C. Lasaga, and D. M. Rye, 1999, Long-term flow/chemistry feedback in a porous medium with heterogeneous permeability: Kinetic controls of dissolution and precipitation: American Journal of Science, v. 299, no. 1, p. 1–68, doi:10.2475/ajs.299.1.1.

Bonnell, L. M., R. H. Lander, and R. Larese, 2013, Prediction of the formation of grain coating chlorite by the in situ alteration of volcanic rock fragments (abs.): 50th Annual Meeting of the Clay Minerals Society, Champaign, Illinois, October 6–10, 2013, p. 20–21, accessed April 4, 2019, http://www.clays.org/50th_annual_meeting_website/docs/CMS-2013-Program.pdf.

Folk, R. L., 1980, The petrology of sedimentary rocks: Austin, Texas, Hemphill, 184 p.

Geloni, C., A. Ortenzi, and A. Consonni, 2015, Reactive transport modeling of compacting siliciclastic sediment diagenesis, in P. J. Armitage, A. R. Butcher, J. M. Churchill, A. E. Csoma, C. Hollis, R. H. Lander, J. E. Omma, and R. H. Worden, eds., Reservoir quality of clastic and carbonate rocks: Analysis, modelling and prediction: Geological Society, London, Special Publications 2015, v. 435, p. 419–439.

Hanor, J. S., 1988, Origin and migration of subsurface sedimentary brines: Tulsa, Oklahoma, SEPM Short Course Notes 21, 248 p.

Hanor, J. S., 1994, Physical and chemical controls on the composition of waters in sedimentary basins: Marine and Petroleum Geology, v. 11, no. 1, p. 31–45, doi:10.1016/0264-8172(94)90007-8.

Hanor, J. S., 1996, Variations in chloride as a driving force in siliciclastic diagenesis, in L. J. Crossey, R. Loucks, and M. W. Totten, eds., Siliciclastic diagenesis and fluid flow: Concepts and applications: SEPM Special Publication 55, p. 4–12.

Harrison, W. J., and L. L. Summa, 1991, Paleohydrology of the Gulf of Mexico Basin: American Journal of Science, v. 291, no. 2, p. 109–176, doi:10.2475/ajs.291.2.109.

Heins, W. A., and S. Kairo, 2007, Predicting sand character with integrated genetic analysis: Boulder, Colorado, Geological Society of America Special Paper 420, p. 345–379.

Hower, J., 1961, Some factors concerning the nature and origin of glauconite: American Mineralogist, v. 46, p. 313–334.

Hutcheon, I., M. Shevalier, and H. J. Abercrombie, 1993, pH buffering by metastable mineral-fluid equilibria and evolution of carbon dioxide fugacity during burial diagenesis: Geochimica et Cosmochimica Acta, v. 57, no. 5, p. 1017–1027, doi:10.1016/0016-7037(93)90037-W.

Land, L. S., G. L. Macpherson, and L. E. Mack, 1988, The geochemistry of saline formation waters, Miocene, offshore Louisiana: Gulf Coast Association of Geological Societies Transactions, v. 38, p. 503–511.

Lander, R. H., and L. M. Bonnell, 2010, A model for fibrous illite nucleation and growth in sandstones: AAPG Bulletin, v. 94, no. 8, p. 1161–1187, doi:10.1306/04211009121.

Lander, R. H., R. E. Larese, and L. M. Bonnell, 2008, Toward more accurate quartz cement models: The importance of euhedral versus noneuhedral growth rates: AAPG Bulletin, v. 92, no. 11, p. 1537–1563, doi:10.1306/07160808037.

Lander, R. H., and O. Walderhaug, 1999, Predicting porosity through simulating sandstone compaction and quartz cementation: AAPG Bulletin, v. 83, no. 3, p. 433–449.

Lee, M-K., and C. M. Bethke, 1994, Groundwater flow, late cementation, and petroleum accumulation in the Permian Lyons Sandstone, Denver Basin: AAPG Bulletin, v. 78, no. 2, p. 217–237.

Le Gallo, Y., O. Bildstein, and E. Brosse, 1998, Coupled reaction-flow modeling of diagenetic changes in reservoir permeability, porosity, and mineral compositions: Journal of Hydrology (Amsterdam), v. 209, no. 1–4, p. 366–388, doi:10.1016/S0022-1694(98)00183-8.

Merino, E., 1975, Diagenesis in tertiary sandstones from Kettleman North Dome. California – II. Interstitial solutions: Distribution of aqueous species at 100°C and chemical relation to the diagenetic mineralogy: Geochimica et Cosmochimica Acta, v. 39, no. 12, p. 1629–1645, doi:10.1016/0016-7037(75)90085-X.

Michard, G., and J. P. Bastide, 1988, Etude geochimique de la nappe du Dogger du Bassin Parisien [in French with English abstract]: Journal of Volcanism and Geothermal Resources, v. 35, no. 1–2, p. 151–163, doi:10.1016/0377-0273(88)90012-1.

Neilson, M. J., and G. F. Brockman, 1977, The error associated with point counting: American Mineralogist, v. 62, p. 1238–1244.

Nesbitt, H. W., 1985, A chemical equilibrium model for the Illinois Basin formation waters: American Journal of Science, v. 285, no. 5, p. 436–458, doi:10.2475/ajs.285.5.436.

Palandri, J. L., and Y. K. Kharaka, 2004, A compilation of rate parameters of water-mineral interaction kinetics for application to geochemical modeling: Menlo Park, California, US Geological Survey Open File Report 2004-1068, 64 p., doi:10.3133/ofr20041068.

Palandri, J. L., and M. H. Reed, 2001, Reconstruction of in situ composition of sedimentary formation waters: Geochimica et Cosmochimica Acta, v. 65, no. 11, p. 1741–1767, doi:10.1016/S0016-7037(01)00555-5.

Paxton, S. T., J. O. Szabo, J. M. Ajdukiewicz, and R. E. Klimentidis, 2002, Construction of an intergranular volume compaction curve for evaluating and predicting compaction loss in rigid-grain sandstone reservoirs: AAPG Bulletin, v. 86, no. 12, p. 2047–2068.

Perez, R. J., and J. R. Boles, 2005, An empirically derived kinetic model for albitization of detrital plagioclase: American Journal of Science, v. 305, no. 4, p. 312–343, doi:10.2475/ajs.305.4.312.

Prieto, M., 2009, Thermodynamics of solid solution-aqueous solution systems, in E. H. Oelkers and J. Schott, eds., Thermodynamics and kinetics of water-rock interaction: Washington, DC, Mineralogical Society of America Geochemical Society, p. 47–85, doi:10.1515/9781501508462-004.

Raffensperger, J. P., 1996, Numerical simulation of sedimentary basin-scale hydrochemical processes, in Corapcioglu, Y. C., ed., Advances in porous media: Amsterdam, Elsevier Science, 440 p., doi:10.1016/S1873-975X(96)80005-9.

Ramseyer, K., J. R. Boles, and P. C. Lichtner, 1992, Mechanism of plagioclase albitization: Journal of Sedimentary Petrology, v. 62, p. 349–356.

Smith, J. T., and S. N. Ehrenberg, 1989, Correlation of carbon dioxide abundance with temperature in clastic hydrocarbon reservoir: Relationship to inorganic chemical equilibrium: Marine and Petroleum Geology, v. 6, no. 2, p. 129–135, doi:10.1016/0264-8172(89)90016-0.

Steefel, C. I., D. J. DePaolo, and P. C. Lichtner, 2005, Reactive transport modeling: An essential tool and a new research approach for the Earth sciences: Earth and Planetary Science Letters, v. 240, no. 3–4, p. 539–558, doi:10.1016/j.epsl.2005.09.017.

Steefel, C. I., and K. T. B. MacQuarrie, 1996, Approaches to modeling of reactive transport in porous media, in P. C. Lichtner, C. I. Steefel, and E. H. Oelkers, eds., Reactive transport in porous media: Washington, DC, Mineralogical Society of America Reviews in Mineralogy 34, p. 83–130, doi:10.1515/9781501509797-005.

Taylor, T. R., M. R. Giles, L. A. Hathon, T. N. Diggs, N. R. Braunsdorf, G. V. Birbiglia, M. G. Kittridge, C. I. Macaulay, and I. S. Espejo, 2010, Sandstone diagenesis and reservoir quality prediction: Models, myths, and reality: AAPG Bulletin, v. 94, no. 8, p. 1093–1132, doi:10.1306/04211009123.

Thyne, G., 2001, A model for diagenetic mass transfer between adjunct sandstone and shale: Marine and Petroleum Geology, v. 18, no. 6, p. 743–755, doi:10.1016/S0264-8172(01)00025-3.

Walderhaug, O., 1994a, Temperatures of quartz cementation in Jurassic sandstones from the Norwegian continental shelf – Evidence from fluid inclusions: Journal of Sedimentary Research, v. 64, no. 2a, p. 311–323.

Walderhaug, O., 1994b, Precipitation rates for quartz cements in sandstones determined by fluid-inclusion microthermometry and temperature-history modeling: Journal of Sedimentary Research, v. 64, no. 2a, p. 324–333, doi:10.2110/jsr.64.324.

Walderhaug, O., 1996, Kinetic modeling of quartz cementation and porosity loss in deeply buried sandstone reservoirs: AAPG Bulletin, v. 80, no. 5, p. 731–745.

Xiao, Y., W. L. Esch, J. E. Welton, R. E. Klimentidis, and P. Rumelhart, 2007, Predicting fault controlled diagenesis and porosity/permeability evolution in Rotliegendes reservoirs, Germany: Insights from reactive transport models: 2007 AAPG Annual Convention, Long Beach, California, April 1–4, 2007, accessed April 4, 2019, http://www.searchanddiscovery.com/abstracts/html/2007/annual/abstracts/lbXiaoYitian.htm.

Xiao, Y., and G. D. Jones, 2006, Reactive transport modeling of carbonate and siliciclastic diagenesis and reservoir quality prediction: 2006 Abu Dhabi International Petroleum Exhibition and Conference, Abu Dhabi, United Arab Emirates, November 5–8, 2006, SPE-101669-MS, 9 p.

Xiao, Y., and G. D. Jones, 2018, Reactive transport modeling and reservoir quality prediction, in Y. Xiao, F. Whitaker, T. Xu, and C. Steefel, eds., Reactive transport modeling: Applications in subsurface energy and environmental problems: Hoboken, New Jersey, Wiley & Sons Inc., p. 157–235, doi:10.1002/9781119060031.ch4.

Xu, T., and K. Pruess, 1998, Coupled modeling of nonisothermal multiphase flow, solute transport and reactive chemistry in porous and fractured media: 1. Model development and validation: Berkeley, California, Lawrence Berkeley National Laboratory Report LBNL-42050, 38 p.

ACKNOWLEDGMENTS

I would like to thank ExxonMobil Upstream Research Company and ExxonMobil Exploration Company for support to conduct this research. The work greatly benefitted from conversations with colleagues Yitian Xiao, Mauro Lo Cascio, David Awwiller, Kathryn Denommee, and Jessica Matthews. The manuscript received valuable and constructive reviews from Rob Lander and an anonymous reviewer, which greatly improved the quality and impact of this paper.

You may also be interested in ...