The Bundelkhand “giant quartz reef” (BGQR) system comprises 20 major quartz reefs which run for tens of km in strike length of average width of 40 m and occurs in spatial intervals of 12–19 km in the Bundelkhand craton, North Central India. The BGQR system is distinct from quartz vein systems originating from crustal scale shearing observed in ancient as well as modern convergent tectonic settings. Fluid inclusions studied in BGQR system are intriguingly diverse although dominated by aqueous fluid which exhibit a broad range of salinity from ~0 to 28.9 wt% NaCl equivalent and temperature of homogenization range of 58 to 385°C. Primary and pseudosecondary aqueous inclusions in assemblages in grain interiors and growth zones vary randomly in their Th—salinity characteristics that preclude identification of discrete fluid events. Aqueous fluid in the BGQR system evolved through mixing of two distinct sources of fluids—a meteoric fluid and a moderate temperature—moderate salinity fluid that was possibly derived from the Bundelkhand granodiorite based on an important clue provided by hydrous mineral bearing fluid inclusions detected by Raman microspectrometry. The results of modeling with PHREEQC indicate that mixing of fluids could be a suitable mechanism in formation of these giant reefs. The available 1-dimensional diffusive transport model for deposition of silica helps in putting constraints on the time span of deposition of silica in the context of the BGQR system. The BGQR system is a possible result of shallow-crustal sources of fluid and silica and could be visualized as a “Paleoproterozoic geothermal system” in a granitic terrane.

“Giant quartz reefs” (Although these quartz bodies have been referred to as “giant veins” by some previous workers, we prefer to use the term “reef” as is done in case of such kilometer scale quartz bodies with/without significant mineralization. Veins, sensu stricto, are embodied in the host rocks and do not impart any geomorphic features on their own unlike the giant quartz “reefs” in the Bundelkhand craton.) are intriguing geologic entities and, as the term suggests, are gigantic veins measuring in tens of kilometers of strike length and tens of meters in width and extend a few kilometers down depth. Such “giant quartz reefs” have been reported from many terranes ranging from Archean to Cenozoic ages such as the Abitibi-Pontiac terrane boundary [1]; Goiás shear zone, Brazil [2]; Bavarian Pfahl, Germany [3]; and South China [4] in a wide spectrum of lithologic ensembles, tectonic settings, and structural locales. These giant quartz reefs are the likely results of mobilization of enormous quantity of silica from suitable source rocks in the earth’s crust and would require fluids of diverse sources and chemistry [1]. Since such processes are primarily controlled by solubility of silica in crustal fluids and operate in widely varying P-T conditions and lithologic environments, they operate in variable time spans in the crust. Precise quantification of these three aspects (source of silica, fluid regime, and time span) is fundamentally challenging. Based on solubility of silica in crustal fluids, deposition of a unit volume of quartz would require two to three orders of magnitude higher volume of fluid [1, 5]. The fate of such huge volume of fluid also needs to be understood in terms of a suitable “sink” that in most cases are assumed to be the country rocks (hydration) hosting the quartz bodies.

A review of literature on formation of quartz (irrespective of form and dimension) in the crust vis-à-vis fluid regime in different domains [612] points to two distinctly different types of models. The first type proposes deeper crustal source of silica and mobilization of the same by fluids generated in the lower crust during metamorphic devolatilization, channelization of the fluid through deep-seated fault/shear zones, and deposition of quartz in form of lodes in response to decreasing pressure-temperature conditions. The seismic pumping vis-à-vis fault-valve model of [9] and the model of “fast-ascending fluid through mobile hydrofracture” suggested by Bons [5] or a variant of it termed as “ballistic” transport of fluid, for formation of quartz bodies of significant dimension, belongs to the first type. The second type visualizes shallow crustal rocks as the source of silica and fluid ascribing prolonged fluid-rock interaction during hydrothermal convection facilitated by enhanced rock porosity/permeability. The first model is often applied by workers on orogenic gold deposits, and the second model applies more appropriately to fossil and live geothermal systems in brittle zones of the crust. Both models are constrained by solubility of silica in crustal fluids to account for the volume proportion of fluid required to deposit the mass of quartz observed in any of the situations mentioned. Kerrich and Feng [1] in their model on origin of giant quartz reefs visualized multiple fluid sources (metamorphic, connate, magmatic, and meteoric) from diverse lithology at variable depths. Occurrence of such giant quartz bodies, irrespective of whether they are on terrane boundaries or in cratonic interiors, poses challenges in establishing the fluid sources and the mechanism of evolution, keeping in mind that fluids of diverse sources have diverse chemistries and variable silica solubilities.

The Bundelkhand “giant” quartz reef (BGQR) system is unique and far outweighs all other occurrences cited above in terms of volume and extent. In the exposed part of the Bundelkhand craton (~29000 km2), there are 20 major reefs [13], each measuring in tens of km in strike length with a nearly uniform interreef spacing of 12 to 19 km that occur all across the craton with an average width of ~40 meters. The major reefs (see Figure 1) and numerous subordinate ones constitute a significant volume fraction of the craton that could only have resulted through extensive mobilization of silica into major planes of weaknesses which are nearly parallelly disposed with a dominant NNE-SSW to NE-SW trend (Figure 1) and subvertical dips.

The present work is an attempt to generate comprehensive data on characteristics of fluid and reconstruct its thermochemical evolution in the process of formation of the BGQR system, based on fluid inclusion microthermometry and Raman spectroscopy. The basic objective of the work is to examine the possibility as to whether such a gigantic system could be a result of a single or multiple episodes of fluid activity or a result of protracted fluid activity over a time period. We have further attempted to discuss the origin of the BGQR system in light of existing models of formation of quartz bodies in the crust on the backdrop of what we decipher on the source of silica and the possible path of fluid evolution and gain insight into the possible temporal span for formation of the BGQR system.

The BGQR system is hosted by the ~2.5 Ga old Bundelkhand granodiorite (BG) [14], which is the dominant phase of the Bundelkhand granitoid complex and constitutes about 80 percent of the exposed area of the Bundelkhand craton (see Figure 1). The TTG gneiss and metamorphosed greenstones exposed mainly around the Bundelkhand Tectonic Zone (BTZ) [13] constitute the pre-BG lithology of Archean age (~3.7 to 2.7 Ga) [13, 15, 16]. Slabunov and Singh [17] reported a range of Pb-Pb zircon dates of ~2.86 to 1.86 Ga from one of the reefs of the BGQR system; the younger dates overlap with deduced ages of minor phases of felsic magmatism in the form of leucogranites and pegmatites [13]. Mafic dykes representing the last phase of magmatic activity cross-cut the Archean units, the granitoid complex, and the BGQR system. Pati [13] quoted an age range of 1.97 to 1.8 Ga for the mafic dykes. The scanty radiometric age data reported so far indicate that the BGQR system was formed between episodes of emplacement of the Bundelkhand granodiorite batholith and intrusion of mafic dykes. The early Proterozoic lithologic ensembles of the Bundelkhand craton were the main source of sediments for the mid-Proterozoic Bijawar and Gwalior basins (see Figure 1) ([18] and references therein). The Bundelkhand craton is fringed by the Vindhyan Supergroup of rocks on the western, southern, and eastern margins. Since the northern margin of the Bundelkhand craton is masked by gangetic alluvium, the juxtaposition of the craton and the Himalayan tectonic zone is obscure. Gokarn et al. [19] ascribed the northward dip of the tectonic fabric of the Bundelkhand craton, north of BTZ, as the effect of Himalayan tectonism. However, there is no clear documentation of any signatures of deformation in rocks of Bundelkhand craton associated with such late tectonic activity.

The occurrence of the BGQR system is elaborated through field features presented in Figures 2(a)–2(f). The salient features of occurrence of the BGQR system are as follows. The major reefs stand out as nearly vertically dipping resistant ridges (as high as 150 meters above ground) resulting in spectacular geomorphic features, identifiable in satellite imageries. Direct contact with the host granitoid as shown in Figure 2(d) is infrequently observed. The width of the major reefs varies between 10 and 70 meters, and many of the major reefs also exhibit apparent variation of width along strike as shown in Figure 2(b). Other than the major reefs, there are numerous subordinate reefs of only in a few hundreds of meters in strike length (a total of ~1500 reported by Pati et al. [13]) that mostly conform to strike of the major reefs. However, a few of them as shown in Figure 2(e) are oblique to the trend of the major reefs. Field features presented in Figure 2 indicate that the reefs of the BGQR system lack prominent foliations and the bulk of it is constituted of white, milky white to light buff-colored quartz. The major reefs have a prominent set of strike parallel joints. Most often, the reefs preserve a stratified nature parallel to the strike of the reef that presumably represent successive growth layers with concomitant deposition of quartz in a gradually increasing fracture apertures. These surfaces and the strike parallel joints facilitated later shearing as depicted in Figures 2(c) and 2(f). Although in most part the reefs are uniformly massive as depicted in Figure 2(h), local scale heterogeneity is imparted by “quartz-in-quartz” type features in the form of multiple veinlets or segregations of white or buff-colored quartz (Figure 2(g)) indicating continued precipitation of quartz in the reefs.

A prominent E-W tectonic zone (the Bundelkhand Tectonic Zone) traverses the craton [13]. The E-W trending Raksa shear zone (RSZ) occurs north of BTZ. However, the temporal relationship of these tectonic zones, especially the BTZ, with the BGQR system is unclear. It may be noted that the BGQR system underwent brittle-ductile shearing (mylonitization) observed in many localities. The EW shear zones significantly offset the strike of the major reefs as observed near Chhatarpur (see Figure 1) area where the reefs attain widths in excess of 50 m. At places sheared surfaces with horizontal striations (lineations) indicating strike-slip movements are observed.

As already stated, the enormous volume of the BGQR system, consisting of reefs trending all across the craton with nearly uniform widths along strike, is unusual and incomparable with occurrences of quartz vein systems globally. The map of BGQR system mostly reconstructed from satellite imageries does not indicate any intricate network of connected reefs and, thus, does not resemble any crustal scale anastomosing shear patterns [20]. Quartz is likely to precipitate from fluid in major structural locales, all the way up from ductile to brittle regimes [21, 22]. The latter favors creation of open space and widening of fracture aperture for growth of veins of larger width. Hydrothermal systems in brittle crust are expected to operate under epithermal conditions and are likely to be longer-lived compared to systems operating at depth under ductile conditions and higher temperature ranges [10, 21]. In the case of these shallow crustal brittle fracture system continuing to greater depths [22], fluids from deeper source can infiltrate to the system. In case of the BGQR system, the down depth extension is largely unknown.

A majority of earlier workers suggested that the BGQR system represents deposition of quartz in prominent (NNE-SSW) and subordinate (at random orientations) fractures as a result of extensional tectonism [17, 2326]. Pati et al. [25] observed sparse chalcophile element (Ag, Sb, Cd, Mo, As, Sn, and Pb) concentrations along with the presence of sulfide minerals in the quartz reefs at restricted localities (in the vicinity of BTZ) and interpreted that a hydrothermal process was responsible for their transport and precipitation. Based on close spatial association with BG, previous workers [24, 25, 2729] ascribed BGQR system to the last phase of granitic hydrothermal activity in the craton.

Representative samples of quartz from the BGQR system were collected in E-W traversing across the major and subordinate reefs covering most parts of the entire craton. Some of the prominent major reefs were sampled along strike to examine intrareef variation in their characteristics. Our emphasis was on larger coverage of samples so as to be representative of the entire BGQR and to obtain a synoptic picture of fluid characteristics. Quartz sampled from the BGQR system was subjected to routine petrography in thin sections (30 μ) under transmitted polarized light using a Leica DM 4500 microscope. The work was primarily aimed at revealing texture of quartz and associated accessory mineral phases and imprints of deformation if any, since same samples were also subjected to fluid inclusion study which is the major component of the work. Fluid inclusion study was carried out on doubly polished wafers (~300 μ) mounted on glass slides using the same microscope. Occurrence of different types of inclusions in the quartz matrix and intergranular and intragranular healed cracks was recorded in the wafers (~25mm×15mm area), and their images were captured at 500x magnification. The selected wafers with workable inclusions were released from the glass slide mounts and fragmented into a number of smaller chips for microthermometry, primarily to accommodate in the sample chamber of the microthermometric apparatus used and also for the ease of locating inclusions during experiments. Heating freezing experiments were carried out using a Linkam THMSG 600 heating-freezing system mounted on a Leica DMLP microscope with 50x long distance working objective. The system was periodically calibrated with synthetic fluid inclusions of pure water and water+CO2 at 0°C and -56.6°C, respectively. Synthetic inclusions of pure water trapped at critical temperature was used to calibrate at 374°C. Heating experiments were carried out only after completion of freezing and recording of all parameters in all inclusion types present in a chip. Reproducibility of phase changes in freezing cycles and heating cycles were ±0.1°C and ±1.0°C, respectively, and inclusions for which reproducibility was poorer were rejected. Heating runs were taken only on those inclusions for which reproducible freezing data could be obtained. The microthermometric parameters recorded for aqueous inclusions are as follows: temperature of first melting (Tfm) (whenever appearance of first melt from a completely frozen inclusion could be clearly observed), temperature of last ice melting (Tm), and the temperature of liquid-vapor homogenization (Th). The temperature of melting of solid CO2 (Tm,CO2) and temperature of liquid-vapor homogenization (Th,CO2) were recorded for the carbonic inclusions. Salinity of aqueous inclusions was calculated from Tm using the equation of Potter et al. [30] and density at Th of the same inclusions were computed using the regression equation of Bodnar [31]. Isochores of the aqueous inclusions were constructed using the formulations of Bodnar and Vityk [32]. Density of carbonic inclusions was computed from equations compiled in [33], and isochores of the same inclusions were constructed using equations from [34] applied to a pure-CO2 condition. All calculations and plotting were done using a Visual Basic Macro enabled MS Excel sheet [35] with implementation of all equations of states in binary and ternary systems corresponding to inclusions studied.

One chip each from samples subjected to microthermometry, with the same population of inclusion types, was used for Raman microspectrometry. A Renishaw InVia Reflex Spectrometer System equipped with a standard confocal microscope was used for Raman spectral analysis. A Renishaw diode-pumped solid-state laser provided 532 nm laser excitation with 5 mW power at the sample. A 2400 grooves/mm grating was used giving a spectral resolution of 0.5 cm-1. Single Raman spectra were obtained using a 5 second integration time with 10 accumulations and a 100x Leica microscope objective, which focused the beam to a spot size of 1 μm. Wavenumber calibration was carried out using an internal silicon standard and was performed as an automated procedure using the Wire version 4.2 software. Two-dimensional Raman mapping was carried out using the StreamLineHR™ mode. The maps were produced by rastering the region of interest in 1 μm steps and recording a Raman spectrum at each step with a 1 second integration time and 1 accumulation. Maps of different minerals were produced by plotting the intensity of definitive Raman bands for each mineral at each step. These were then displayed as colored maps using the Wire software.

On the backdrop of the picture of fluid characteristics obtained from fluid inclusion studies, geochemical modelling was attempted using PHREEQC Interactive 3.6 (USGS software repository) software package. The 1-dimensional model of dissolution and diffusive transport of silica in crustal fluid was used with incorporation of advection of fluid to visualize the possible timeframe of formation of the BGQR system. Calculations of different parameters from the analytical solution to the dissolution-diffusion-advection equation were performed, and their graphs were generated using MATLAB R2017a.

4.1. Petrography of the BGQR System

The most commonly observed petrographic features in the BGQR system is depicted in Figures 3(a)–3(d). The dominance of quartz is well observable while other minerals such as K-feldspar and mica/sericite constituting around 5-10% of volume of the reefs. Occasionally, sericitized K-feldspar and plagioclase, rare chloritized biotite, epidote, and zircon occur as isolated patches in quartzose matrix that are remnants of slivers of hydrothermally altered host granitoid. There is remarkable variation in grain size from very fine to very coarse-grained nature, sometimes blocky (Figures 3(a) and 3(b)). A typical mosaic of coarse quartz grains is shown in Figure 3(a) with grains mutually impinging the growth of each other indicating deposition of quartz in open space. Figure 3(b) shows vary coarse grain aggregates with feeble extent of cataclastic deformation (peripheral granulation) with the quartz megacrysts showing undulose extinction.

4.2. Fluid Inclusion Types

Fluid inclusion study was focused on quartz as described in the previous section with negligible imprints of deformation as shown in Figures 3(c) and 3(d). As revealed from petrography (Figure 3), the BGQR system preserves pristine fluid deposited textures. The chosen fluid inclusion population in quartz in BGQR system is likely to preserve their pristine compositional characteristics in the host quartz under such situation as demonstrated by Tarantola et al. [36]. In the present study, microthermometric runs were taken keeping the possible effect of later deformation in mind, restricting to inclusions with moderate to small (between 5 to 15 μ in maximum dimension) sizes and without any visible distortion of shape to be the least affected in any post entrapment deformation. Thus, we make a reasonable assumption on the reliability of our microthermometric data representing the fluid as regards to the thermal and salinity characteristics.

Fluid inclusion features such as occurrence of inclusions in three-dimensional network and distribution in healed fractures as trail bound inclusions and each with the tentatively identified compositional type were recorded. We have essentially discriminated fluid inclusions based on occurrence in random three-dimensional framework, spatially away from healed fractures (from now onwards, simply “random three-dimension”) and those on transgranular and intragranular healed cracks in host quartz in order to understand the different pulses of fluid activity covering the whole spectrum of fluid inclusions in different mode of occurrences. We have focused our attention to the majority occurring randomly in 3 dimensions in the host quartz that satisfy the criteria of primary inclusions [37], yet, we do not presume them to have been entrapped at the same instant of time and at the same P-T-X condition of the fluid. We have not been able to definitively “finely” discriminate events of fluid inclusion entrapment [38] in the absence of cathodoluminescence (CL) images of the quartz grains, hence, cannot describe them as discrete FIAs. Wherever growth zones are easily demarcated in a quartz grain under the microscope, discrete FIAs were considered and examined their distribution in the temperature-salinity space that we have achieved after the microthermometry. However, as acknowledged by Chi and Lu [39] and Chi et al. [38], it is commonly not possible and easy to establish FIAs in all cases. Besides, the criteria for discriminating FIAs do not mention about the spatial restrictions in particular, and it is not certain whether FIAs could be correlated across fields of view in a sample as well as across samples. In a long-lived hydrothermal system, it is more likely to be a case of continued precipitation/recrystallization of host minerals (especially quartz) in which case there could be multiple FIAs which eventually have to be taken as a whole to understand the path of fluid evolution. We have adopted the approach as suggested by [40] (Figure 3 of the author) for mode of growth of hydrothermal quartz resulting from protracted fluid activity with concomitant fracturing. Fluid entrapped at any stage of growth (as primary inclusions) is also likely to be entrapped as “trail bound” inclusions in healed cracks in the immediately preceding stage and all could be a part of the continuum of fluid activity. Therefore, in our case, inclusions present in healed cracks that generally assigned as secondary could be part of the same continuum. Since our approach is to present a scenario of the entire spectrum of fluid activity, we have focused mostly on inclusions distributed “randomly in three dimensions.” Even if all such inclusions might not have been trapped at the same time, rather at different times in recrystallized domains (but undeformed), we are still in a position to get the scenario of changing fluid characteristics in terms of temperature and salinity distribution of data. Further, in such protracted fluid activity in the brittle domain in the crust, minor variations in chemistry at different stages and local scale heterogeneity due to mixing of fluids are expected [41]. Therefore, occurrence of inclusions of different chemistry in the same field of view may not always indicate entrapment at different times and may be part of the same assemblage. We have considered a small representative population of inclusions in healed cracks (trail bound) for a comparison with the major population occurring “randomly in three dimensions.”

A few representative of the veins (as shown in Figure 2(g)) possibly resulting from the protracted hydrothermal activity were also included for fluid inclusion studies to check for any variation in them. Thus, in the absence of CL imaging, a clear picture of broad temporal events can be achieved with respect to the fluid inclusions in the giant quartz reefs as main event to that of subordinate quartz vein as the latest stage of its formation. A total of thirty-two wafers (28 from the quartz reefs and 4 samples representing quartz veins in the reefs) were selected based on abundance of fluid inclusions. Based on the contents (phases) of fluid inclusions at room temperature determined by their response to freezing-heating cycles, a total of five inclusion types are recognized. They are type I: aqueous; type II: pure carbonic; type III: aqueous carbonic; type IV: solid phase bearing aqueous polyphase inclusions; and type V: unusual. A synoptic description of the various fluid inclusions types is furnished as follows.

4.2.1. Type I: Aqueous Biphase Inclusions

Type I inclusions are the most dominant in all samples of the BGQR system. There is a scarcity of workable type I inclusions in the fine-grained matrix of quartz. In the case of later veins, in which distinct growth zoning is observed infrequently, there is also a dearth of fluid inclusions. The samples with observable growth zones display alternate milky and transparent quartz where mostly irregular nature of inclusions are seen in the milky band and good workable inclusions in the transparent band (Figure 4(l)). The type I inclusions occur in their typical biphase (L+V) state at room temperature as shown in Figure 4(a). Their size varies from <2 μm to 30 μm. Elliptical to oval shapes are dominantly seen followed by shapes like rectangular, square, rhombic, and polygonal to irregular ones. Negative crystal shape of some of these inclusions is also a notable feature. The irregularity of shape of these inclusions seems to be a function of size. Type I inclusions often occur in monophase liquid only condition which may be ascribed to their entrapment at very low temperature-pressure and simply lack a vapor bubble due to metastability [42]. Association of type I inclusions with any of the other types is observed in clusters as well as in transgranular and intragranular trails.

4.2.2. Type II: Pure Carbonic Inclusions

These inclusions occur in biphase (L+V) state at room temperature with “pseudo-Brownian” movement of a vapor bubble that helps in their identification. Size ranges from <2 μm to 8 μm with oval to elliptical shapes (Figure 4(b)). In all instances, they occur in association with type I inclusions (Figure 4(j)). A few of these inclusions occur in monophase (liquid only) state and could only be confirmed with appearance of vapor on cooling and response to freezing at very low temperatures.

4.2.3. Type III: Aqueous Carbonic Inclusions

This type is the least abundant in the BGQR system and occurs as both biphase [Aq (L)+Carb(L)] and tri-phase conditions [Aq(L)+Carb(L)+Carb(V)] (Figure 4(c)). Their size varies from <5 μm to 15 μm: average 8-10 μm with elliptical and oval shapes. They occur either isolated or in association with type I and/or type II inclusions in clusters (Figure 4(k)).

4.2.4. Type IV: Aqueous Polyphase Inclusions

Aqueous inclusions are at times observed to contain solid phases with tabular, square, prismatic, rhombic, or oval (Figures 4(d)–4(f) and 4(l)) habits. Although they are quite frequent in occurrence in the studied samples, the solid phases in them are quite varied in their color and refractive index. Such inclusions also occur without a vapor phase in the inclusion cavity. These inclusions vary in size from <5 μm to 20 μm. In some instances, it is seen that nearly half of the inclusion cavity is occupied by a bright phase as shown in Figure 4(f). The solid phases in the type IV inclusions do not respond to heating.

4.2.5. Type V: Unusual Type Inclusions

This type of inclusions is essentially aqueous in nature with presence of a spherical vapor bubble and typically contains a single or multiple tiny objects that execute continuous random movement within the aqueous liquid in the inclusion cavity as shown in Figures 4(g)–4(i)(movies of the behavior of these inclusions at room temperature have been uploaded on https://www.youtube.com/channel/UCo1PDrZtrdMLS9Lqv3rHeDw). These objects are of variable sizes (filling up most of the inclusion), and their vigorous movement imparts an “effervescence” type appearance with the entities colliding among themselves. The fast random movement would normally indicate that these entities have lower density than water, but not vapor. In case they are liquid hydrocarbon, they can only be pentane or heavier, and in that case, the density contrast with liquid water would be too insignificant to cause such random movements. The other possibility could be that the tiny entities are solid mineral particles under some special situations since any mineral should be denser than the aqueous liquid and would normally be static in the liquid area as in all polyphase inclusions. Careful observation reveals that they are quite frequent in many of the samples of the BGQR system. They occur exclusively as clusters or occasionally in association with type I inclusions. Although random movement of these entities is quite common, often they also occur without executing any motion. On some instances, these inclusions appear to contain a glassy phase along with other solid phases with or without a vapor. We categorize inclusions with a suspected glassy phase as “suspected melt inclusions” and have discussed them in Section 6. All the aforementioned compositional fluid types are essentially observed in the “random three-dimensional network” within the quartz matrix of quartz reefs which constitute variable compositional fluid inclusion assemblages as tabulated in Table 1. It can be noted here that the type II, type III, and type V inclusions are absent in the trail bound inclusions as well as in the subordinate quartz veins. Moreover, the presence of type II and type III CO2-bearing inclusions in the matrix quartz is also rarely encountered in that indicate negligible proportion of carbonic fluid, at least by volume during the formation and evolution of the giant quartz reefs. If at all any discrimination in light of assemblage is possible in or case, we could include the rare type IV and more frequent type V inclusions into one assemblage and the aqueous biphase (type I) inclusions into the other, with subassemblages of coexisting type I and type II inclusions formed under locally varying environmental parameters. A further analysis is presented in the following section with the help of Figure 5. The fluid inclusion population in all samples of BGQR is dominated by type I inclusions (see Table 2). Figure 4 depicts all modes of occurrences observed in the BGQR system. Sample-wise occurrence of inclusions of all types is summarized in Table 2.

4.3. Microthermometry

The following discussion is based on microthermometric data on 854 fluid inclusions from the bulk of quartz making up the reefs and 85 fluid inclusions from subordinate quartz veins in these reefs. It is to be noted that the following discussion is essentially focused on the fluid inclusions in the random three-dimensional networks targeting within the quartz grains as mentioned in the previous section. They are labelled as matrix quartz in the subsequent data plots. These are compared with the fluid inclusions of other mode of occurrence, viz., trail bound in the matrix quartz and subordinate quartz veins in the quartz reef. Although five different types of fluid inclusions are inferred in the petrography section, the current discussion on the microthermometric data is restricted to only type I and type II inclusions. This is because of the incomplete retrieved dataset from the type III, type IV, and type V inclusions in terms of salinity information and complete homogenization temperature. The type III inclusions decrepitated before total homogenization due to buildup of internal pressures. Solid phases in the type IV inclusions also did not show any indication of dissolution up to their decrepitation. In type V inclusions, no significant phase changes were observed except that the rate of movement of the particles slowed down with the lowering of temperature below 0.0°C. They become static in the completely frozen aqueous liquid in the inclusion. They get back to their previous state after melting of ice and further rise in temperature. In such inclusions, Tm varies in a narrow range from -0.2 to 0.0°C having salinity value 0.0 to 0.35 wt% NaCl eq. Such low values of salinity of the aqueous liquid indicate that the moving objects, if solid, are captive phases. The partial homogenization of the inclusions is achieved with the disappearance of vapor bubble in a temperature range from 191 to 267°C corresponding to densities of 0.76 to 0.88 g/cc. There is no indication of disappearance of these vibrating entities with further heating. The inclusions eventually decrepitate on heating in excess of 400°C. Thus, the overall interpretation of the microthermometry work is dominantly based on the type I inclusions with some type II inclusions. A precise description is as follows:

A broad range in Tm and Th value is evident when quartz grains with visible growth zonings are observed with the help of petrological microscope. An attempt has been made to represent them schematically (Figure 5: not to the scale) to depict the data distribution for Tm and Th, in view of their variable modes of occurrences as discussed in the petrographic section following [43]. Figure 5(a) is a simplified illustration of a portion of a grain with growth zones which represent a situation in the matrix quartz domain of major quartz event of formation of the giant reefs. A thorough insight on the nontrail bound inclusions in the core part represents a situation where near pure water (Tm ~0.0 to -4.2°C) with a moderate to higher temperature values (Th ~162 to 330°C) are evident, whereas the 1st growth zone (subsequent to the core) depicts a relatively more saline fluid type (Tm ~-13°C) with a lower temperature signature (Th ~125°C). Though such low temperature (Th ~128 to 137°C) values are observed in the 2nd growth zone, the Tm is found to be of comparatively higher value (i.e., Tm ~-8.8 to -9.5°C) to that of first.

Additionally, it comprises moderately higher temperature values (Th ~157 to 196°C) with pure water (Tm ~0.0 to -0.1°C) signature as is seen in the core part of the quartz. Moreover, the fluid inclusions in the third growth zone preserve even more saline type (Tm ~ -24°C) with very low-temperature value (Th ~100°C) which occur in association to what observed in the first zone and some from the core. When these fluid inclusions which are essentially part of the 3D network (or sole occurrences) are compared to that of some intragranular trail bound inclusions in the same domain, they furnish qualitatively either similar event or no clear to distinct event of fluid activity which can be deciphered unambiguously from them to establish any chronological order. In other words, the fluids entrapped in the core as well as the growth zones do not represent any distinct event of fluid activities rather more akin to interplay of fluids of differing salinity and temperatures resulting in a broad salinity and Th ranges of aqueous inclusions. In contrast to Figure 5(a), Figure 5(b) represents a part of the subordinate quartz veins in the quartz reef showing the mesoscopic zonings. The growth zones present in both sides of the chip represent the growth direction of the quartz grains. Hence, the fluid inclusions present in the central portion of the schematic diagram (Figure 5(b)) represent the latest phase of fluid activity. Although they largely represent a moderate salinity value (with Tm ~-13 to -15°C) and Th of 130 to 150°C, a critical examination reveals the following: (i) a comparatively lower Tm value of about -15°C is associated with Th of <150°C, whereas Tm ~-13°C can be correlated to >150°C; (ii) there is a broad variation in Tm of about -13 to -17°C within a narrow Th value of ~150°C in different FIAs; (iii) additionally, very few inclusions with a near zero value of salinity are also observed with a broad Th value, i.e., from 167 to 298°C in different FIAs, whereas a higher Tm of about -21°C and a low Th value of 130°C is also seen to occur as alone in the central portion of the chip; (iv) the two of pseudosecondary fluid inclusion assemblages signify almost similar Th values (~170°C) with a very contrasting salinity, viz., Tm ~-0.6 and -15.4°C, while the secondary trail represents Tm and Th of about -20°C and 150°C, respectively; and finally, (v) a secondary trail (youngest) cross-cutting the whole growth zones as well as final projecting quartz grains represents the Tm value of -13°C to -16°C and Th value of 140°C to 150°C. These values are easily correlated to the values seen in the variable growth zones present in both the sides. Even these growth zones which predate the secondary trail are also associated with some very lower values of Tm ~-18°C to -21°C and also very low Th ~84°C to 120°C. Thus, such a broad variation as well as similarity in the variety of mode of occurrences with the observed salinity and Th cannot be a result of reequilibration event rather interaction of at least two different compositional fluid types, viz., low salinity (~0-2 wt% NaCl eq.) and moderately high salinity (~15-20 wt% NaCl eq.), is evident.

The temperature of first melting (Tfm) could be recorded only in about 7% of the type I inclusions, which furnished a broad range from -20.8 to as low as -61.4°C and is irrespective of salinity value. This indicates the presence of multiple cations other than Na+ though the number of fluid inclusions to show such behavior is very less. Besides this, brown coloration on complete freezing due to formation of antarcticite [37] is observed in many such inclusions which furnish Tfm values that indicate presence of appreciable Ca in the fluid. The recorded final ice melting temperature in the presence of vapor (Tm) ranges from ~0.0 to -30.0°C. Salinities calculated from Tm values show a broad range from a lower value—almost pure water to a higher value of 28.9 wt% NaCl equivalents (Figure 6(a)). The clustering is mostly around low salinity values. More than half of the type I inclusions give salinity values within 2 wt% NaCl eq. A similar type of data range is also inferred for the fluid inclusions when variable domains are considered as per our classification (matrix quartz, trail bound, and quartz veins in the reef). Although the distribution of data points is bimodal in case of the subordinate quartz veins, the dominance of low saline inclusions is well observable around 0-2 wt% NaCl eq. with subsequent higher values afterwards in all the cases. The temperature of homogenization (Th) of type I inclusions ranges from 72.5°C to 385°C as shown in Figure 6(b) with rare values below 100°C or above 400°C. The Th data of type I inclusions mostly cluster around 150 to 250°C, about 11% of it giving values exceeding 250°C. Microthermometric data are summarized sample wise in Table 2, as per locations shown in Figure 1. All type I inclusions homogenize to liquid phase without a single instance of vapor phase homogenization indicating that the aqueous fluid never underwent boiling. This indicates that entrapment was from a homogeneous liquid (irrespective any event of fluid activity) at all stages of formation of the BGQR system. Density calculated using the regression equation of [31] ranges from 0.582 to 1.191 g/cc (Figure 6(c)). Inclusions from samples of quartz veinlets have a dominance of high density inclusions.

Type II inclusions furnish a near constant value of Tm,CO2 values at -56.6°C indicating that the carbonic phase is devoid of any other volatile species such as CH4 or N2. The temperature of homogenization (Th,CO2) of these inclusions ranges from 7.8°C to 30.9°C (Figure 6(d)), and the mode of homogenization is always to liquid phase furnishing densities in the restricted range of 0.53 to 0.88 g/cc (Figure 6(e)).

4.4. Raman Microspectrometry

Raman microspectrometry was attempted on selected representative fluid inclusions of different types. It may be noted that considering the extensive nature of BGQR and the diversity of inclusion types, a comprehensive and exhaustive Raman spectroscopic characterization has not been possible.

The type I aqueous biphase inclusions were scanned routinely for detection of molecular species in both the liquid and vapor phases. Raman analysis of the liquid phase showed the characteristic broad OH stretching bands between 2800 cm-1 and 3800 cm-1 (Figure 7). As reported in [4446], the shape of these bands varies with salinity. In agreement with the microthermometry data, the majority of the Raman spectra showed two broad OH stretching peaks at approximately 3230 cm-1 and 3400 cm-1 (Figure 7(a)) indicating a low salinity aqueous phase. A smaller proportion of inclusions showed a narrower OH stretching band at approximately 3450 cm-1 (Figure 7(b)) indicating fluids of higher salinity which broadly agree with salinity data of type I inclusions. Due to the relatively large amount of noise in the Raman spectra, no attempt was made to determine the salinity of the fluid using the methods cited above. The Raman spectrum of one type IV inclusion also showed the presence of an additional OH stretching band at 3552 cm-1 (Figure 7(c)). This has been tentatively assigned to the presence of a solid inclusion believed to be chamosite. Raman analysis of type II and III inclusions confirmed the presence of CO2 as shown in Figure 8(a) No other gases could be detected by Raman spectroscopy. Raman bands at 282 and 1087 cm-1 also confirmed the presence of calcite in some type IV inclusions (Figure 8(b)).

Type V inclusions are the most intriguing ones, barely reported in literature. Raman microspectrometry was attempted on these inclusions. It has not been easy to identify the fast-moving objects under the beam; however, mapping has been able to reveal their characteristics to some extent. No hydrocarbon species could be detected in these inclusions. Raman maps of selected type V inclusions are shown in Figure 9. Figure 9(a) shows an image of a type V inclusion containing many moving particles. A Raman spectral map of the same inclusion is shown in Figure 9(b). Raman spectra containing the main OH band of clinochlore (at 3571 cm-1: Figure 9(b)) are shown by the red pixels. This indicates that the majority of the solid phases in this inclusion are clinochlore. Another solid phase with a Raman band at 3367 cm-1 (see Figure 9(g)) is shown by the green pixels in Figure 9(b) and is identified as epidote. Some of the other solids remain unidentified, but the Raman spectra did not indicate the presence of any hydrocarbons in this inclusion. According to Bodnar (R. J. Bodnar, personal communication), the tiny mineral particles being hydrous silicates (clay) absorb the IR component of the microscope illumination and execute random movement as they are thermally agitated. The presence of such tiny mineral particles either as captive or daughter phases makes an interesting case as regards to the source of fluid for these BGQR as discussed in a later section.

A suspected melt inclusion with a small vapor bubble is shown in Figure 9(c). The Raman map in Figure 9(d) shows that the glassy phase has bands at 285, 331, and 512 cm-1 (Figure 9(c)) and is anorthite (green pixels in Figure 9(d)). The red pixels in Figure 9(d) correspond to the main Raman band of CO2 at 1385 cm-1 and confirm that the vapor contains CO2. Another suspected melt inclusion is shown in Figure 9(e) and contains a brownish solid and a clear, partially recrystallized phase. The Raman spectra of this inclusion shown in Figure 9(f) contain bands at 478 and 507 cm-1 (Figure 9(g)) indicating the presence of albite (green pixels in Figure 9(f)). A weak Raman band is also observed at 700 cm-1 (Figure 9(g)) and is assigned as mica (red pixels in Figure 9(f)). Figure 9(g) embodies the representative spectra of the solid phases present in these inclusions as described in Figures 9(a)–9(f).

A limited idea about the possible pressure-temperature variation of fluid in the BGQR system can be obtained from Figure 10(a) where isochores of a few coexisting type I and type II inclusions are plotted. We assume that they represent coeval inclusions and attempt to estimate the trapping temperature and pressure conditions. Isochores of type I and type II inclusions were constructed in P-T space using available formulation for NaCl-H2O and pure CO2 fluids. The pairs of isochores intersect in P-T space furnishing temperature and pressures of 240°C and 1409 bars, 244°C and 678 bars, 260°C and 857 bars, 263°C and 635 bars, and 364°C and 1820 bars (Figure 10(a)). This method of thermobarometry proposed initially by Roedder and Bodnar [47] assumes of simultaneous entrapment of the CO2 and aqueous liquids in the host mineral. Barring two of the high-pressure values which may not satisfy the condition of immiscibility in the H2O-CO2±NaCl system [34], rest of the values, although scanty, could be considered. These ranges indicate fluctuating conditions of pressure, possibly between purely lithostatic to intermediate situations corresponding to depth ranges of 3 to 4 km. A hydrostatic pressure condition at such depths (as low as 300 to 400 bars) would indicate opening of the major fractures to the surface. One important thing that seem to emerge is that we get a crude idea about the pressure correction that our Th data would require. If we take average pressures of 600–800 bars from the three acceptable intersection points, the Th values of aqueous inclusions would need a correction of 25–40 degrees; thus, peaks in our Th histograms would shift to 200–220°C.

A distinct presence of a very low saline fluid (0-2 wt% NaCl; temperature varying from less than 100° C to 250° C) on Figure 10(b) indicates that a meteoric fluid endmember is conspicuous. Higher temperatures in the same low salinity range (~350°C) would possibly indicate but not definitive clue to a metamorphic fluid of deeper source due to the CO2-poor nature. On the other hand, the high-temperature signature could also be a result of transient heating of the same meteoric fluid. Near uniform distribution of points across a range of salinity (from near zero to 30 wt% NaCl eq.) within a narrow variation of temperature (~150 to 250°C) is also quite striking. Such variation of salinity can be ascribed to mixing of fluids as is observed in many cases (c.f. [48]) and also demonstrated through numerical modelling [49]. Although fluid mixing most likely operated in the formation of the BGQR system, a definitive identity of endmember components can only be speculated since fluid inclusion characteristics represent the regime of deposition due to mixing of the endmembers. Therefore, mixing of a meteoric water (likely to be colder than the observed minimum temperature) with another fluid endmember of moderate salinity and temperature (likely higher than the observed limits) seems to be the path of evolution of fluid in the BGQR system. The moderate temperature-salinity fluid endmember could represent either a residual fluid from the host Bundelkhand Granodiorite (BG) that evolved to such low temperatures through internal evolution or evolved to moderate temperature salinity through cooling and dilution by meteoric fluid and was present as a pervasive fluid before the brittle deformation episode. This could be the reason of different salinity and temperature signatures of fluid types observed along the core of the host grains as well as in the growth zonings (Figure 5). Even a fluid of similar temperature value witnessed a broad salinity value or the vice versa. Such internal evolution of granitoid derived fluid or rather presence of lower temperature-salinity fluid within an evolving granitic pluton with involvement of meteoric fluid has long since been discussed in the literature [48, 50, 51]. Although low-temperature and moderate salinity fluids in granitoids causing alterations are usually ascribed to incursion of meteoric water [51], in some instances, a single magmatic fluid evolving to low-temperature salinity state through fluid-rock interaction has also been demonstrated [50].

Rout et al. [52] observe striking resemblance in fluid characteristics between the BGQR system and that in BG. This provides the first clue to link the origin of the BGQR system with BG in terms of fluid regime. The second important clue is provided by occurrence of inclusions with tiny particles of hydrous minerals (identified as clinochlore, epidote, and unknown hydrosilicate phases as of now) in the BGQR system (described as type V inclusions and discussed at length) as well as in the BG domain, where they occur in greater frequency [53]. Such types of unusual inclusions are observed in shallow crustal epithermal systems (R J Bodnar, per. Comm.) in which there is fluid flow from country rock to fracture zones in response to drop in pressure. We visualize intermittent advective flow to fracture zones which were the locales of deposition of quartz, and the fluid carried these mineral particles in suspension. It may be noted that Wilkinson et al. [54] observed anhydrous minerals within aqueous inclusions and identified a “silicothermal” fluid as product of felsic magmatism. In our case, based on the hydrous nature of captive phases, we do not visualize the fluid to have evolved through such a mechanism since we do not have a clue that such captive phase bearing fluid was of any higher temperature. We visualize a moderate temperature aqueous fluid that evolved internally in BG and advected to the fracture zones in a pressure gradient, carried tiny particles of the hydrous minerals (those formed during deuteric alteration of primary minerals in BG), and deposited quartz.

The temperature-salinity characteristics indicated from microthermometric data and depicted through the fluid evolution diagram do not compare well with metamorphogenic fluid typified by the well-studied quartz veins in orogenic gold deposits as regards to the insignificant CO2 component, that too, without any CH4 and also the overall pressure-temperature-salinity ranges [55, 56]. The BGQR system therefore appears to be dominated by meteoric fluid akin to geothermal systems in brittle deformation zones in the crust. The fracture zones provided the locales for mixing of fluids derived partly from the host BG (present in fracture networks or the pervasive fluid in intergranular spaces) and partly from deep circulating meteoric fluid. The advection of fluid from the country rock is likely to have taken place along a pressure gradient developed due to fracturing in an extensional phase, depending on the magnitude of pressure difference [7, 10]. Slabunov and Singh [17] also suggested E-W extension of the Bundelkhand craton as the cause of extensive fracturing and formation of the giant quartz reefs. The major reefs likely evolved through a concomitant increase in the fracture aperture and deposition of silica over a period of time as long as the extensional phase would have lasted. Such successive growth layers provided easy planes of shear as illustrated in Figure 2(f), observed commonly in these reefs. Increase in fracture aperture has alternatively been ascribed to pressure (mechanical work) exerted by growing crystals from a supersaturated solution on the fracture wall [5759]. Our fluid inclusion data do not point to multiple episodes of fluid activity since we do not observe any disparity in fluid characteristics between the major reefs, minor quartz veinlets that likely formed as a result of protracted fluid activity. Temperature-salinity characteristics of fluid in quartz associated with sparse sulfide mineralization in BG [28] also conform to the scenario presented here. Therefore, minor metal enrichment in the fluid is a more likely result of local chemical control rather than different episodes of fluid activity. Although it is not possible to work out a proper temporal framework for such a process of formation of the BGQR, a qualitative idea can be obtained through numerical modeling as discussed in the subsequent section.

The main objective of the modelling exercise is to examine the suggested path of fluid evolution inferred from fluid inclusion studies. It is imperative to have the insight into the mineralogical control of solubility of silica in a shallow crustal condition. Further, it is also important to see whether the proposed mixing of fluids could be an efficient process in deposition of quartz. We also visualize a scenario of advection of fluid derived from the host BG and mixing on fracture zones that formed the “giant reefs” of the BGQR system. Therefore, we make an attempt to extend the 1-dimensional model of dissolution-diffusion transport proposed by Wangen and Munz [8], to see whether the process could facilitate formation of the quartz reefs.

6.1. Modeling with PHREEQC

We attempted simple models for examining the solubility of silica in the fluid equilibrated with a granite mineralogy at 300°C and used a starting composition of the fluid with concentrations of ions (Na+, K+, Ca2+, and Cl-) in a geothermal fluid and equilibrating with a modal granite mineralogy (primary as well as secondary phases such as chlorite and kaolinite) using the LLNL database of PHREEQC. This fluid represents the moderate-temperature endmember as discussed before. In the absence of precise elemental/ionic analyses of inclusion fluids, the temperature and composition of the fluids taken for the modelling exercise is only a crude approximation based on our results of microthermometry. However, with the SIT database, a dramatic increase in solubility of Si is observed with anorthite as an equilibrium phase (values as high as 2.2 mole/kg when anorthite is increased to 5 moles). Our results are expected to have some uncertainties due to extrapolating to temperatures higher than the range of SIT database. We have not taken kinetics into consideration in both cases. It is observed that the SIT thermodynamic database results higher solubility of Si due to higher concentrations of H3SiO4-, H2SiO42-, and Ca (H3SiO4)+ species, which are not accommodated in the LLNL database. In order to proceed a bit further in the modeling exercise, an attempt was made to examine the effect of mixing of fluids for obtaining a semiquantitative evaluation of the formation of the BGQR system. The PHREEQC input file was designed for mixing of the two fluids—the granitoid derived fluid at 300°C with the acquired chemistry on interaction with a granite (with primary and secondary mineral phases). The meteoric fluid is taken at a temperature of 50°C at pH of 7 and ideally free of any dissolved electrolyte. Figure 11(a) shows results of such exercise where the mixing was carried out at 9 : 1, 8 : 2, 7 : 3, 6 : 4, and 5 : 5 proportions. The corresponding variation of pH of the mixed fluid is shown in Figure 11(b). Our attempt of modelling of fluid mixing is in line with Walter et al. [41] who also presented varying saturation indices of minerals with mixing proportions of two fluids with constraints from chemistry of inclusion fluids and directly sampled fluids. The scenario of increasing saturation index of quartz with increasing proportion of meteoric fluid remains the same with LLNL database, differing only in the concentration of silica.

As can be seen from the plot, saturation index of silica increases with increasing fraction of the cold meteoric water with the temperature of the resultant fluid decreases to about 180°C, and correspondingly, the molality of silica decreases in the fluid mixture. Although speculative, the results of this modeling exercise clearly indicate the efficiency of mixing with cold and dilute meteoric fluid in giving rise to deposition of silica (Figure 11).

The high values of solubility obtained with SIT database, if real, would far exceed the experimentally measured solubility even at very high pressures [60]. This would imply that we may not have to invoke a deeper source of silica for formation of quartz reefs, and additionally, the volume requirement for mobility of silica could be greatly reduced. Dove and Nix [61] deduced a higher solubility of Si in the fluid in the presence of cations such as Na+, K+ and Li+, and Ca++ and Ba++, than in pure water based on kinetic considerations. Initial solutions with appreciable Li+ and variable concentrations of Na and K were also considered in the model with a granite mineralogy. The results obtained do not indicate any appreciable increase in solubility of Si. Since clay minerals are believed to be enhancing silica solubility [62], different fractions of clay minerals (Na and Ca montmorillonite, kaolinite, and also chlorite) of variable mass fractions were included in the equilibrium assemblage. Solubility only up to 0.25 moles/kg is observed in the presence of kaolinite. In case the maximum solubility of silica in a fluid that interacted with an anorthite-present mineralogy can go to 2.2 moles per kg (about 130 g/kg of water) and the fluid on mixing with meteoric fluid at 50°C undergoes a drop in silica solubility to ~1.2 mole/kg (~70 g/kg of water), it would result in deposition of 60 g of silica that translates to about 22 cc of silica per 1000 cc of water. In other words, for each cc of silica deposited, it would require about 45 cc of water, as compared to 240 cc of water in case of fluid of deeper origin [5].

At this stage, a commentary on the PHREEQC modelling may be as follows. The high solubility obtained by using the SIT database is too high compared to values reported at 300°C taking crystalline silica only as the equilibrium phase [63]. Mysen [64] observed a reduced solubility of silica in deeper crustal fluid in presence of Ca-rich basic rocks corresponding to deeper crust. In that context, high solubility obtained in the present case at 300°C in the presence of anorthite needs to be confirmed. In summary, although the result of the exercise attempted with SIT database looks promising for formation of BGQR type system with silica derived from shallow felsic crust, a more rigorous analysis is needed for a more definitive quantification.

6.2. Dissolution-Diffusion-Advection Model

The mixing model discussed in the foregoing section indicates that formation of quartz on fracture zones in shallow crustal environment through fluid mixing is feasible. Such a process implies physical movement of fluid (advective transport of silica) from the source that needs to be examined from a space-time perspective. Wangen and Munz [8] proposed a model of 1-dimensional dissolution and diffusive transport of silica that helps in visualizing the space-time framework, other than the important fact that the model does not require a huge proportion of fluid: quartz. Although the model was proposed to explain quartz veins of smaller dimensions in metamorphic rocks, it possibly could be examined in the present context. The model visualizes enhanced solubility of fluid with respect to quartz in the fluid in country rock matrix due to pressure solution. The dissolution and diffusive transport process was stated in a nondimensional form and analytically solved for a steady-state condition by bringing in the Damkohler number (Da), a nondimensional parameter representing the relative dominance of dissolution-diffusion and precipitation (see [8] for detailed derivation of the model). Wangen and Munz [8] proposed a rate of cementation of quartz in the fracture as a function of the Damkohler number and fundamental kinetic parameters and degree of supersaturation of the fluid with respect to fluid in the fracture. They considered the characteristic time and lengths of diffusion and dissolution to be same and used the Damkohler number (Da) to represent the relative dominance of diffusion-dissolution and precipitation and showed the normalized concentrations at different normalized distances at fixed degree of supersaturation at different values of Damkohler number (varying in four orders of magnitude: 0.01, 0.1, 1, and 10 showing increasing dominance of diffusion). As expected, the values are high at low Da. Further, they also showed the relationship of rate of cementation (dw/dt) as a function of temperature and Da and showed nearly six orders of magnitude variation in rate of cementation as temperature varies between 300°C to 50°C and at least 4 orders of magnitude variation in the rate with more than 3 orders of variation in Da. Our results of PHREEQC modelling qualitatively imply that with mixing of fluids, the rate of cementation increases with increasing proportion of the cooler component. As stated before, there are minor issues that were not given much importance by [8] such as the equilibrium constant of the reaction of dissolution of SiO2 did not take the dependence of density of water into account that has been incorporated in our calculation using available formulations [31, 63]. Such corrections are not expected to make significant change to results in terms of steady-state concentration profile since the formulation takes normalized concentration. However, it is expected to make some difference while calculating the rate of cementation as discussed later. Since the system under consideration involves flow of fluid to the fracture zone in response to change in pressure, advective transport of silica in the fluid could not be ignored, and therefore, we expanded the formulation (equation (1) of [8] with incorporation of an advection term (see equation (A.1), Appendix, where the terms are as in [8]). The third term on the L.H.S. of the equation represents advection (v—velocity in m/s). Boundary conditions imposed are same as in the original formulation, that is, a Dirichlet boundary condition at x= at any time t (m=m2) and a Neuman type boundary at x=0 at any time t (equation (A.2), Appendix). Nondimensional distance and time were introduced as the original formulation and equation (A.1) was reformulated in its nondimensional form (equation (A.3), Appendix) where τ is the normalized time, c is the normalized concentration, and χ is the normalized distance. The Dirichlet boundary condition is represented by c=1 at χ=, and the Neuman boundary (equation (A.4), Appendix) is represented, where Pe is the Peclet number (=vl0/D). The system was further reduced to a steady-state condition and solved by standard technique [65]. Figure 12(a) shows the normalized concentration profile generated using equation (A.8) (Appendix) taking Pe=0 to check the validity of our analytical solution. It may be noted that the analytical solution of [8] eliminated the degree of supersaturation. We retained the term since the parameter controls the concentration gradient for the diffusive transport. However, the concentration profile remains the same with values of normalized concentration at any normalized distance varying only in the fourth place from decimal with increase of the degree of supersaturation. The profiles of normalized concentration with higher contribution of advection (higher Peclet number) are shown in Figure 12(b). It is observed that advection does not have a positive effect on the normalized concentration at any distance. Saishu et al. [11] showed higher rate of deposition of quartz with dominance of advection that is expected to be reflected in higher normalized concentration at any distance compared to a “diffusion only” situation. However, the present advection-diffusion transport model fails to indicate that. Figure 12(c) is a contour plot of normalized concentration in a Damkohler number (Da) versus Peclet number (Pe) plot which is a better depiction of the combined effect of diffusion and advection. Incorporation of advection in the equation for the cementation rate of Wangen and Munz [8], was done with Peclet number resulting in equation (A.10) (the appendix). Figure 12(d) depicts the scenario of the cementation rate versus temperature for different values of Da without any advection contribution. Here, all parameters (degree of supersaturation, effective surface area, molar volume of quartz, and activation energy) are same as that of Wangen and Munz [8]. However, we have taken a different value of the prefactor kinetics of quartz (2.02×101 rather than 6.3606 mole m-2S-1) from [66]. The figure shows time span of formation of quartz vein of a unit half width (1 m) with diffusion only as the transport mechanism (Pe=0). The figure indicates that for deposition of a 1 m half width quartz vein it would take 2.043 million years at a temperature of 300°C corresponding to the lowest Da of 0.01 and is longer if higher values of Da are taken. Since it is a one-dimensional approximation, it is presumed that the time of formation holds good to a quartz vein of any length and thus applies to situations pertaining to BGQR. A reef of 20 m half width would take about 40 million years to form. It would be important to note that here, the effect of degree of supersaturation is quite prominent on the rate of cementation. At a value of the degree of supersaturation S=0.1 (one order of magnitude higher than the value taken by [8], the time span of deposition of a half width of 1 m is reduced to 0.2 million years (see Figure 12(f)). Further, if we take a depth of 3 km in the crust and consider the pressure difference of what exactly is the difference between purely lithostatic and hydrostatic (about 45 M Pa), the cementation rate is increased by three times and the time span by reduces to one third, that means, at 3 km depth, a quartz vein of 1 m half width will take 0.074 million years to form, implying that a 20 m half width (average of 40 m width for any major reef of the BGQR system) will take about 1.5 million years.

As stated already, incorporation of advection does not improve the situation of cementation rate as shown in Figure 12(e) using equation (A.10) (Appendix) where contours of cementation rate are shown on a Da versus Pe plot at a fixed temperature. It is interesting to note that there is hardly any increase in the cementation rate with introduction of advection, at a fixed Da. As can be seen in the figure, at any Da, cementation rate decreases with increasing Peclet number. However, lowering the Damkohler number requires higher Peclet number to maintain a particular cementation rate. It may be noted that Saishu et al. [11] visualized dominance of advection atΔP>10 MPa and considered pressure difference up to 160 MPa. They did not explicitly mention the values of Peclet number in their calculation although the same parameter in hydrothermal flow can be as high as 106 ([10], who use molecular Pe as the parameter that is applicable to such situation). The system possibly needs to be formulated with different other parameters as discussed in [67] who critically evaluated such nondimensional parameters (especially the Damkohler number) and opined that they are of limited applicability. In a diffusion only model, one order of magnitude higher in degree of supersaturation results in an order of magnitude increase in the cementation rate and thus would result in reduction in time of formation of a unit half width of quartz vein (keeping the effective surface area same) which is of significance to the BGQR system. Wangen and Munz [8] used the ratio of the characteristic length of dissolution-diffusion (ld) and characteristic length of precipitation (lp) (see their equation (12)). The equation indicates that the value of Damkohler number is controlled by the value of concentrations of silica in the country rock and the equilibrium concentration of fluid in the fracture and the porosity since other parameters (diffusion coefficient, kd, kp) are fixed at a particular temperature and fluid chemistry (mainly pH). For example, at 25°C, with known values of these three parameters (the kinetic parameters taken from Palandri and Kharaka, 2004), a value of 0.005 (0.5%) for porosity, and value of specific surface area of fracture walls as unity, the value of Da is ~0.08 and there is no scope of orders of magnitude variation of the value since the only parameter to influence Damkohler number is porosity. Wangen and Munz [8] took a porosity value of 0.5%. An order of magnitude decrease in the value will result in a Da value close to 1, but the permeability would be too low to be realistic. Moreover, a Da value of the order of 0.01 is already low indicating negligible control of diffusion on deposition of quartz in the fracture fluid and thus introduction of the advection term in the present form of the governing equation remains ineffective. Therefore, a precise quantification of the process needs refinement of the model for silica transport in fluid.

Before going on to draw conclusions on the origin of the BGQR system, it would be worthwhile to examine the possibility of deeper source of fluid and silica for the BGQR system. The “mobile hydrofracture” model of [5] proposed a fast ascent (1 m/s) of fluid from deep source at pressure in excess of 5 kb and temperature of about 600°C, so as to inhibit any loss of silica during its upward migration and deposit quartz at a shallow depth only in response to drop in pressure and temperature. Bons [5] visualized the fluid to get arrested at a depth of about 6 km in 3 hours, yet lose heat and decrease in temperature to 250°C. If the fluid follows an isothermal path (as shown in Figure 1 of [5]), it is highly unlikely that the movement of the fluid could be as fast. A slow upward migration of the fluid would continuously precipitate quartz on its upward migration with decrease in both pressure and temperature. On the other hand, if the fluid actually moves as fast as 1 m/s, it would be following an adiabatic path. A crude estimate indicates that from a pressure of 5 kilo bars to a pressure of 2 kilo bars, the temperature would drop from 500°C to 460°C and the solubility of silica would drop to less than half and would deposit quartz. Such high-temperature deposition of quartz should be manifested in fluid inclusion homogenization temperatures. The range of fluid inclusion temperature obtained from BGQR would not agree with such a scenario even applying a pressure correction as discussed before. In the light of discussion by Spear (1993) (see Figure 19-19 of [21]) it would be difficult to visualize fracturing and creation of so much of open space for transport of fluid and silica in voluminous quantity at depths beyond 5 km. It may be noted that Staude et al. [68], in their model of extension-driven dewatering of middle crust, assumed prevalence of hydrostatic conditions to depths down to 13 km in their model. However, they clearly stated that their model involves large uncertainties and numerous poorly constrained parameters in geological, physical, and chemical processes. We prefer the depth range of prevalence of hydrostatic conditions to be 5 km and find agreement in the paucity of quartz reefs of such large dimension (as the BGQR system) in orogenic gold provinces where seismic pumping of fluid from deeper source is inferred. Further, and most importantly, the model also does not explain convincingly the fate of the large volume of fluid after deposition of quartz. The model of [5] and resultant style and geometry of deposition of quartz (Figure 11 of [5]) does not agree with the observed mode of occurrence of the reefs in the BGQR system; hence, it does not seem tenable to explain the formation of the BGQR system.

According to Saishu et al. [11], the time taken to seal a crack is greatly reduced with advection of fluid from the rock matrix to these cracks and advection dominates with increase in ΔP. They used ΔP values of 0 to 25 MPa and observe that at small values up to 0.7 MPa, advection is low (due to low ΔP) and the process is diffusion controlled and the rate of deposition of silica in cracks is quite low. Advection is dominant at higher value of ΔP, and at 25 MPa, the time taken for sealing a crack with aperture of 56 μm and length of 7.4 cm is about 56 years with a porosity of 3% in the source rock. If we consider the half width of 20 m for each reef of the BGQR system and similar magnitude of pressure difference to have operated, then each reef would have taken about 40 million years (the strike length remaining unknown) to form that possibly would compare well with results of numerical modelling involving transport of silica in the system by diffusion only. The process depicted in [11] quite explains the phenomenon in a subduction zone (convergent tectonic regime) and quartz veins in multiple numbers is visualized due to opening of new fractures in subsequent seismic rupturing rather than fluids being channeled through preexisting sealed fractures (single pass flow through a fracture rather than multiple pass). It would be difficult to extrapolate such a situation to the BGQR to explain the fracture zones occupied by the reefs to have resulted from discrete seismic rupturing events over a period of time. Williams [69] discussed the ineffectiveness of coseismic boiling of fluid in deposition of quartz and quoted a rate of ~1 μm of deposition of quartz per kilo year that translates to 1 mm/Ma which is far too slow compared to values calculated using the 1-dimensional diffusive transport (low Da situation) and precipitation model. Therefore, processes taking place in a convergent setting involving seismic rupturing and fault-valve mechanism cannot be similar to formation of giant quartz reefs in a cratonic interior and may involve different mechanism of fluid-rock interaction. According to [70], the cementation rate is a much more complex phenomenon and is quite variable in the case of formation of quartz veins in meter scale under hydrothermal conditions of concomitant fracturing and growth of quartz. Within the limitation of the data generated in the present work, we can only qualitatively discuss some alternative models of origin of BGQR.

Although we find the model of [8] plausible to explain our case, its applicability in toto is fraught with an inherent problem. Any diffusion modelling exercise is essentially targeted to profile the concentration in space and time using the standard Fickian diffusion equation in form of a parabolic partial differential equation. Such models are usually applied to transport of any chemical species in a certain medium (solid, liquid, or gas) in space (considered in one dimension) with reference to a point source and concentrations at all other points at the initial time (t=0) is taken as zero. In the case of transport of silica from country rock and deposition in open spaces, the scenario is likely to be different. If we consider an initial representative elementary volume (width x) representing the half width of the fracture space (Figure 13(a)), an initial batch of fluid occupies the fracture and is interfaced with the fluid in the country rock in immediate contact. Deposition on the fracture surface would set up the concentration gradient for subsequent diffusive transport of silica to the fluid in the fracture. If the country rock is discretized to small volume elements as shown, the diffusion of silica in the fluid should be taking place across infinitesimally small distances across each volume element and finally at the interface of two fluids at the fracture. In that case, the characteristic time and lengths of diffusion and dissolution and the Damkohler numbers do not seem to be quite applicable. The process will be mostly controlled by the kinetics of quartz dissolution and precipitation and will be controlled by the one which is slower.

We visualize that major reefs of the BGQR system occupy near-vertical extensional fractures. If we make a crude estimate of the percentage of space that need to be created for these fractures, it turns out to be about 0.1% or 0.001 (if we take a length of 400 km for the width and the total length of fracture aperture to be 400 meters). Such low degree of extension is very unlikely to be a result of any lithospheric scale extensional tectonism. Hippertt and Massucatto [2] suggested a feasible mechanism of creation of open space as a result of loss of volume on shear planes being compensated by the extension gashes to be filled by fluid and in turn quartz veins. They described km scale quartz veins in such tension gashes between crustal scale mylonite zones in host granite. The situation presented in [2] is that of a Pan-African reactivation (strike-slip tectonics) of structures in a Proterozoic granite giving rise to extensional gashes and formation of the quartz veins. In case of the BG (~2530 Ma), there is no documentation of such parallel disposition of intense and wide mylonitic zones. In case such zones exist, the major reefs of the BGQR system should be at high angles to such zones for which there is no clear documentation. Lastly, the extension gash veins studied by [2] are of smaller dimension compared to BGQR and would indicate much smaller fraction of volume created for deposition of quartz. In the absence of a detailed documentation of deformation in BG, such a possibility remains open. As mentioned earlier, such fracture zones provided the locales for later shearing as observed in the craton and also manifested in the BGQR. A more rigorous modelling is needed to substantiate this as one of the tentative propositions on the origin of the BGQR.

This leaves us with the alternate proposition of a “shallow crustal source” for the fluid and silica and the most feasible source for both is the host BG. We could visualize a prolonged episode of extensional fracturing [13, 17], of whatever possible cause that created zones of reduced pressure to have ingress of fluid from the BG surrounding each major fracture. Such a fluid infilling of the fracture zone triggered deposition of quartz and the process sustained with diffusive-advective transport of silica to the fracture (Figure 13(b)). Bons et al. [71] suggested that the downward flow of (meteoric) fluid could be brought about by desiccation of fluid by hydration of rock. Such a process could possibly have operated that gave rise to pervasive alteration of minerals in BG. However, in the absence of features of exhumation, we prefer a lateral flow of fluid in response to fracturing and reduction of pressure. In such a scenario, continuous advection of fluid from the surrounding could be possible only if there was concomitant increase in the fracture aperture since the fracture is unlikely to have been as wide at the beginning as the present width of the reefs. Although not fully supported by results of numerical experiments involving advection, flow of supersaturated fluid from the surrounding BG could be one of the mechanisms for the inferred continued deposition of quartz along with contribution from diffusive transport. Fluid inclusions, especially the particulate hydrous silicate bearing inclusions in quartz in the host BG and comparable characteristics of fluid in both domains [53], lend credence to such a proposition. Our proposition on the origin of the BGQR system compares favorably with the evolution of the quartz lode in the Bavarian Pfahl system [3, 72, 73] and the fault-controlled quartz reef in the South China bloc [4]. Both these systems involve shallow crustal derivation of fluid and silica and fluid-rock interaction as the dominant process although the two differ significantly in terms of the time span of hydrothermal activity. A comparison of fluid characteristics with continental live geothermal systems may be worthwhile. The two live geothermal systems associated with granitic terranes are the Broadlands-Ohaaki geothermal system in Northland, New Zealand, and the Kakkonda geothermal system in NE Japan. The BGQR system broadly resembles the Broadlands-Ohaaki system in terms of the temperature and salinity ranges of deep fluid, albeit absence of boiling of aqueous fluid in BGQR. The Kakkonda geothermal system is mostly a liquid dominated system having a deep (>300-350°C; ~3-4 km) and a shallow (230-260°C; <1500 meter) reservoir [74, 75]. Muramatsu and Komatsu [76] reported precipitation of hydrothermal quartz from the upward flow of fluid in the productive fractures at shallow depths. K-Ar dating from hydrothermal K-feldspar in equilibrium with the quartz also suggest the evolution of the Kakkonda geothermal system after the emplacement of the Kakkonda granite which might have acted as the main heat source for the geothermal system [77]. Doi et al. [78] described the Kakkonda geothermal system to be strongly controlled by neotectonic fractures. The Kakkonda geothermal system is a CO2-poor system comparable to the BGQR system. The broad salinity and Th values are very much comparable to each other which varies from 222-380°C and 0-29 wt% NaCl+CaCl2 in case of Kakkonda geothermal system. The Kakkonda system has been visualized as a result of mixing between meteoric fluids with deep fluids. The similar interpretation is also made from the current work in the BGQR system within limitation of the work elements. Even the coeval type I and type II inclusions in the BGQR suggest a similar depth range of about 3-4 km for their formation with a prevailing temperature of ~240 to 260°C as is for the deep reservoir of the Kakkonda geothermal system. Moreover, like the Bundelkhand granitic complex, the Kakkonda granite also represents a composite of quartz porphyrite, tonalite, granodiorite, quartz diorite, and granite. Based on these broad comparisons and fluid characteristics depicted and combined with results of modelling, we can visualize the BGQR system as an ancient “geothermal system” that operated in the shallow crust.

The BGQR is a result of fluid activity in the shallow crust, within a granitic terrane, with fluid mainly linked to the late-stage fluid that evolved internally within the evolving granitic pluton or a late pervasive fluid there. Fluid inclusion temperature-salinity characteristics in the BGQR nearly rule out deeper fluid sources. Occurrence of inclusions with suspended hydrous silicate phases (chlorite and epidote) raises the possibility of advective flow of fluid from the host granite. Fluid mixing appears to have been the dominant mechanism of deposition of quartz as a result of decreased solubility of silica in response to emplacement of fluid in locales of reduced pressure and also cooling brought about by mixing with cooler meteoric fluid. Extension of the 1-dimensional dissolution and diffusive transport model [8] with incorporation of advective transport of silica to the fracture zones provides important insight on the possible time span of formation of the BGQR system. Although a crude estimate of the time span of formation of the giant quartz reefs could be made from the present modelling exercise, a more accurate picture needs a proper evaluation of contributions from dissolution, advection, diffusion, and precipitation. The model also needs refinement with a more realistic mechanism of the transport-deposition process. The very disposition of the BGQR in the craton indicates that they are unlikely to arise from crustal scale shearing process such as those operating ancient orogenic belts (terrane boundaries). Processes operating in present-day active convergent tectonic settings cannot possibly be extrapolated to explain the BGQR system. Rather, BGQR is more likely a Paleoproterozoic analog of live geothermal systems in a granitic terrane with involvement of meteoric fluid. Silica for the BGQR system is likely sourced from the host BG, possibly in more than one way—from the late stage evolved silica-rich fluid that migrated into the fracture zones and the meteoric fluid, also carrying silica in its downward/lateral flow through fracture network. Both these fluid components deposited quartz on major and subordinate fractures. Although there are not many present-day analogs of such an ancient system, the geothermal field in Northland, New Zealand [6], possibly comes close. A direct comparison with any live geothermal system may not be conclusive due to the fat that the older systems like BGQR are remnants after significant unroofing. Second, there is no such granitic terrane of batholithic dimension to compare with BG. A more rigorous analysis is needed for a proper understanding of this ancient system, in terms of the nature of fluid flow and source of heat. More importantly, absence of significant mineralization in BGQR also needs to be accounted for. Lastly, it would also be worthwhile to state that understanding the exact mechanism of creation of such major fractures zones with near uniform spacing in a Paleoproterozoic granitic complex (BG) would certainly give important insight to the evolving tectonism of the lithosphere, not only in the Bundelkhand craton but also in the global context.

Appendix

A. Advection-Diffusion Model for BGQR System

The governing equation on incorporation of the velocity (advection) is given by
(A.1)ϕmtxϕDmx+vϕmx=kd1mm2,
where the notations are as in Wangen and Munz [8]. ϕ represents porosity of the rock, and m, m1, and m2 represent concentrations (in mol/m3) concentration of silica in the fluid at any point, concentration of silica in the fluid in the fracture, and concentration of silica at infinite distance from the fracture, respectively. The third term introduces advective transport with v representing horizontal velocity boundary condition of m=m2 at x= holds. The constant flux boundary condition at x=0 can be stated as
(A.2)ϕDmxϕvm=kpamm11.
Using the same scheme as Wangen and Munz [8], nondimensional concentration (c), time (τ), and distance (χ) were derived with characteristic time of diffusion (l0) and dissolution (ld) and characteristic time of diffusion (t0), dissolution (td), and precipitation (tp). The 1-D advection-diffusion equation was formulated as
(A.3)cτ2cχ2+Pecχ=t0td1c,
where Pe is the Peclet number given as vl0/D (v—velocity m/s). Dirichlet boundary condition is as before, that is, c=1 at χ=, and the Neuman boundary condition can be stated as
(A.4)cχ=Pe+t0tpc+PeSatχ=0,
where S=m2m11. The characteristic time of dissolution and precipitation are as in Wangen and Munz [8].
Equation (A.3) was restated for a steady-state condition as
(A.5)2cχ2Pecχ=Nd1c,
with the Dirichlet boundary condition imposed as before and the Neuman boundary at χ=0 as
(A.6)cχ=Pe+Npc+PeS.
By using a separation of variable method (Kumar, 2007), the analytical solution to equation (A.5) with the imposed boundary conditions is derived as
(A.7)c=1+Pe+Np+PeSPe+Pe2+4/2Pe+NpexpPe+Pe2+42.
Following the same logic as in Wangen and Munz [8], with characteristic time of diffusion and dissolution being same, Np could be replaced by the Damkohler number for the system Da, and equation (A.7) is restated as
(A.8)c=1+Pe+Da+PeSPe+Pe2+4/2Pe+DaexpPe+Pe2+42.

Equation (A.8) is used to calculate the normalized concentration of silica at varying normalized distance with variable Pe and Da.

We have used equation (22) of Wangen and Munz [8] for the rate of cementation of quartz on the fracture and reformulated their equation (23) as given below incorporating Da and Pe.
(A.9)mfm11=S1+Pe+Da+PeSPe+Pe2+4/2Pe+Da.
Finally, we derived the equation for the cementation rate as stated below:
(A.10)dwdt=2kpavqS1+Pe+Da+PeSPe+Pe2+4/2Pe+Da.

There is no additional data on this manuscript.

The work presented here is a part of Duttanjali Rout’s Doctoral Thesis.

The authors declare that there is no conflict of interest in regard to publication of this paper.

DR acknowledges the host institute for financial assistance in the form of Research Assistantship during the tenure of her doctoral program and the host department for infrastructural support. TM acknowledges that this work has been made possible by access to the Raman Lab at the ACT node of the Australian National Fabrication Facility (ANFF – ACT). MKP acknowledges support from the Department of Science & Technology, Government of India and Ministry of Earth Science, Government of India, for past funding for the equipment support for fluid inclusion work and the host institute for extending financial support for the computational facilities. We thank Pritam Nasipuri, Editor of this Special Issue of the Journal, for inviting us to contribute and also for his valuable editorial support.

Exclusive Licensee GeoScienceWorld. Distributed under a Creative Commons Attribution License (CC BY 4.0).