Evidence of Hydrocarbon Generation and Overpressure Development in an Unconventional Reservoir Using Fluid Inclusion and Stable Isotope Analysis From the Early Triassic, Western Canadian Sedimentary Basin

Deep burial of sedimentary basins results in the development of complex diagenetic environments influenced by pressure, temperature, and metasomatic chemical processes. Fracture systems resulting from deep tectonic-related burial can provide archives of physio-chemical characteristics during burial helping unravel diagenetic events such as hydrocarbon migration and paleobarometry. The Early Triassic Montney Formation in the Western Canadian Sedimentary Basin is a highly productive unconventional hydrocarbon reservoir that has undergone multiple phases of tectonic-related burial and uplift resulting in the formation of a series of calcite-filled fracture systems. These fracture systems occur as vertical to sub-vertical fractures, brecciated zones, and horizontal bedding-plane parallel fractures that are rich in co-occurring, but not co-genetic aqueous and petroleum fluid inclusion assemblages. Fluid inclusion microthermometry, Raman spectroscopy, and stable isotope analysis of these fracture systems and host rock reveals paleobarometric and temperature conditions during fracture formation. Vertical fractures formed at temperatures exceeding 142°C during peak burial associated with the Laramide orogeny ∼50 Ma. Similarities in modeled oxygen isotope values of calcite parent fluids and pore water implicate locally sourced carbonate in fracture calcite. Therefore, low permeability and closed system-like conditions were prevalent throughout initial fracture formation and cementation. Petrographic analysis of brecciated and horizontal fractures show evidence of hydrocarbon generation and migration into fracture-filling calcite. Modeling of petroleum inclusion paleobarometry indicates entrapment pressures approaching or even exceeding lithostatic pressure consistent with the development of overpressure associated with the thermal maturation of organic matter following peak burial. Combined use of aqueous and petroleum fluid inclusions in this deeply buried sedimentary system offers a powerful tool for better understanding diagenetic fluid flow, the timing of hydrocarbon migration/maturation, and helps constrain the pressure-temperature history important for characterizing economically important geologic formations.


INTRODUCTION
Unravelling the diagenetic history within deeply buried sedimentary systems is an essential component to understanding how processes such as tectonic activity, fluid migration, and hydrocarbon maturation evolved in the basin. These processes can result in the formation of fractures providing records of the physio-chemical conditions during their formation. Fractures in sedimentary rocks form as a result of stress and strain, related to the response of a rock undergoing lithostatic, tectonic, thermal stress, or high-pressure fluids (e.g., Hubbert and Willis, 1957;Cosgrove, 1995;National Research Council, 1996;Lai et al., 2021). Due to the intrinsic relationship between fractures and economically important petroleum and mineral systems, assessing their composition and characterizing parent fluids within these systems is essential for understanding their role in subsurface fluid flow.
Fluid inclusions can play a vital role as they provide pressuretemperature and geochemical snapshots of these vein forming fluids, including primary vein forming fluids, but also subsequent fluids that migrate through these fractures (Roedder, 1984). Parent fluids captured within inclusions provide direct evidence of fluid composition, while microthermometric analysis can provide evidence of pressure-temperature conditions during entrapment (e.g. Roedder, 1984;Goldstein, 2001;Dubessy et al., 2001;Goldstein et al., 2003;Burruss, 2003;Bakker, 2018). This evidence, coupled with burial history models, can also help elucidate the timing of entrapment (Gao et al., 2017;Gasparrini et al., 2021). Therefore, detailed analysis of fluid inclusion assemblages can unravel complex diagenetic histories within structurally complex systems. In addition, petroleum fluid inclusions provide valuable archives of the properties, distribution, and evolution of hydrocarbons in the subsurface and when used in conjunction with aqueous fluid inclusions can provide pressure-temperature history of their host reservoir (Pironon and Bourdet, 2008;Volk and George, 2019).
The lower Triassic Montney Formation is one of the most productive unconventional reservoirs in Canada (National Energy Board, 2013). As a consequence, it has been cored and studied extensively, providing a large collection of material for analysis. However, owing to its complex diagenetic history (Liseroudi et al., 2020(Liseroudi et al., , 2021Vaisblat et al., 2021), regional inconsistencies, and variable tectonic history (Wozniakowska et al., 2021), more work is needed to better characterize this economically important resource. Unconventional systems associated with low-permeability conditions exhibit closed system-like behaviour (Cesar et al., 2020) and often can be overpressured. Fluid inclusions provide a history of formation pressure and temperature and are therefore instrumental in understanding past diagenetic phases in unconventional petroleum reservoirs (Pironon and Bourdet, 2008;Gao et al., 2017). To date only one study has evaluated fluid inclusions in fractures within the Montney Formation (Gasparrini et al., 2021), however extensive work has been conducted on fluid inclusions within the Devonian and Mississippian strata of the Western Canadian Sedimentary Basin (e.g. Al-Aasm and Clarke, 2004;Al-Aasm et al., 2019 and references within; Davies and Smith, 2006;Yang et al., 2001). Devonian strata exhibit widespread evidence of hydrothermal alteration via the occurrence of saddle dolomite and high temperature alteration of organic matter (Davies and Smith, 2006), distinguishing characteristics that are not observed in the Montney Formation in British Columbia. Therefore, more work is needed to better constrain this diagenetically active system. Lastly, deciphering diagenetic history aids in distinguishing diagenetic from primary geochemical signatures, which is essential in the evaluation of paleoenvironmental versus postdepositional alteration geochemistry including primary versus migrated hydrocarbons (e.g., Ardakani et al., 2020b).
In this study, we use a combination of fluid inclusion microthermometry, stable isotope geochemistry, and Raman spectroscopy to evaluate the fluid sources and pressuretemperature conditions in which calcite-filled fractures formed in the Montney Formation. Modeling of the aqueous and petroleum inclusions reveals the physio-chemical characteristics of the system including the sources of parent fluids and development of overpressure. Finally we assess the implications for the evolution of hydrocarbon maturation and migration in the Montney Formation.

GEOLOGICAL SETTING
Sample material for this study was collected from a core of the Montney Formation within the Western Canadian Sedimentary Basin (WCSB) (Figure 1). Development of the abundant mineralized fractures within this core may be associated with one of several post-depositional events between the mid-Triassic and Eocene such that a brief overview of the history of the WCSB is provided for context.
The WCSB is a westward-thickening sedimentary wedge that is up to 6 km thick adjacent to the fold-and-thrust belt of the present-day Canadian Rocky Mountains in the west; and which thins to an erosional edge in the east on the Canadian Shield (Price, 1994;Wright et al., 1994;Tufano and Pietras, 2017). The fold-and-thrust belt of the Rocky Mountains of the Canadian Cordillera was shaped over four tectonic pulses from Late Jurassic to Late Cretaceous, each corresponding to a depositional pattern in the nearby foreland basin (Pană and van der Pluijm, 2014). Paleogeographically, the Montney Formation of the WCSB was deposited on a clastic ramp along the northwestern margin of the Pangea Supercontinent (Golonka, 2007) in the Peace River Embayment (PRE) ) in a collisional retro-foreland basin (Rohais et al., 2018).
The Montney Formation unconformably overlies the Late Permian Belloy Formation and is unconformably overlain by organic-rich shales of the Middle Triassic Doig Formation, or the Lower Jurassic Gordondale Member of the Fernie Formation in British Columbia and Alberta, respectively (Davies et al., 2018;Zonneveld and Moslow, 2018). In the distal western part of the basin, the Montney Formation entails offshore transition and offshore organic-rich siltstones as well as turbidite deposits where the Montney Formation is overlain by the Lower Doig organic-rich phosphatic zone Davies et al., 2018).
The Montney Formation is divided into three informal submembers: the Lower Montney, Middle Montney, and Upper Montney Members separated by the Dienerian-Smithian and Smithian-Spathian stage boundaries (Davies, 1997;Davies et al., 2018;Zonneveld and Moslow, 2018). The Lower Montney Formation occurs throughout the WCSB and primarily consists of fine-to medium-grained, laminated bituminous dolomitic siltstone. The Middle Montney Formation represents a thick succession of bituminous dolomitic siltstone and minor fine-grained sandstone with a low abundance of bioturbation. The Upper Montney consists of a thick succession of fine-to coarsegrained, bituminous, dolomitic siltstone, and it is completely removed due to erosion throughout much of the Alberta portion of the basin (Zonneveld and Moslow, 2018).
Deposition of the Montney Formation was significantly influenced by paleostructural elements in the Peace River area including: the underlying Upper Devonian Leduc reefs, the Hines Creek-Fort St. John Graben, and northwest-southeast-trending faults such as the Dunvegan and Tangent faults . Normal faults, including the Carboniferous to Permian Dawson Creek Graben Complex and precursor Precambrian basement normal and strikeslip faults, are abundant in the Peace River Embayment (PRE) Davies et al., 1997;Hope et al., 1999;Berger et al., 2008).
It has been seismically and structurally demonstrated that extensional normal faulting was episodically reactivated from the Devonian until the Cretaceous (Hope et al., 1999;Mei, 2009), particularly by loading of compressional thrust faults associated with the Jurassic Colombian and Cretaceous Laramide orogenies (O'Connell et al., 1990).
Diagenetic evolution and compositional modifications in the mineral assemblages of the Montney Formation occurred significantly from Early Triassic (early stage corresponding to shallow burial) to Late Cretaceous (Late stage corresponding to deep burial; Liseroudi et al., 2021;Vaisblat et al., 2017Vaisblat et al., , 2021 and prior to maximum burial (80-60 Ma, Ness, 2001;Davies, et al., 1997;Ducros et al., 2017;Rohais et al., 2018). The Montney Formation has undergone uplift and erosion from the end of Paleocene to the present time (Ness, 2001;Ducros et al., 2017;Rohais et al., 2018). The Montney Formation entered the hydrocarbon generation window at approximately 90 Ma (Ness, 2001). The timing of maximum burial of the Montney Formation in British Columbia is estimated to be between 50-60 Ma (during the Laramide Orogeny) with a maximum burial temperature of approximately 160-170°C, adequate for thermogenic gas generation (Ness, 2001). This is evident from the increasing Montney thermal maturity trend to the west and southwest (Wood and Sanei, 2016;Gibbs and Rakhit, 2019;Euzen et al., 2021).

Sample Material
Sixteen calcite-filled fracture samples were collected from a single core near the western margin of the Montney Formation in WCSB (Figure 1). Four doubly-polished thick section slides (100 µm thick) were prepared from the largest calcite fractures ( Figure 2). This included two samples from the Lower Montney, one sample from the Middle Montney and one sample from the boundary of the Upper Montney and overlying Doig Formation (Figure 1). Vertical fractures were selected as opposed to horizontal fibrous "beef" calcite due to the small size of fluid inclusions in horizontal fracture calcite. Samples were prepared using a liquid cooled rotary saw to prevent modification of fluid inclusions.

Petrology and Fluid Inclusion Microthermometry
Fluid inclusion petrology followed the fluid inclusion assemblage (FIA) concept (Goldstein and Reynolds, 1994) to divide fluid inclusions into groups that were trapped contemporaneously. Detailed petrology was conducted on a Zeiss Axio Scope A1 Microscope fitted with a Linkam THMSG 600 Temperature control stage to heat to over 200°C and cool below -100°C. Petrographic images were captured with a Zeiss AxioCam ICc 5 camera using both white and ultraviolet light. Petrographic work was carried out at the Organic Petrology Laboratory, Geological Survey of Canada, Calgary. Calibration of this system was conducted using synthetic H 2 O and CO 2 fluid inclusion standards for a precision of ± 1°C at 375°C and ± 0.2°C at -56.6°C. Homogenization temperatures were measured first to prevent stretching or decrepitation of fluid inclusions due to ice formation. Salinity was calculated from final ice melting temperatures (Tm) using the salinity-ice melt relationship of Bodnar (1993): Salinity (wt.% NaCl) 1.78 · ( − Tm) − 0.0442 · Tm 2 + 0.000557 · (Tm3). (1) Qualitative determination of the electrolyte composition in aqueous fluid inclusions was assessed by measuring the eutectic temperature (Te) when fluid inclusions devitrified following heating of frozen inclusion (Bodnar, 2003).
Isochore and bubble point curve calculations of aqueous fluid inclusions were modeled via the AqSo_NaCl program (Bakker, 2018) using measured homogenization and last ice melting temperatures. Petroleum phase envelopes and vapor fraction curves were calculated using the PVT simulation package PVTsim (Calsep Co.) using a recombined hydrocarbon composition derived from locally produced fluids (data from field operators). Isochore calculations for petroleum fluid inclusions were modeled using PVTsim software using petroleum fluid inclusion homogenization temperatures, petroleum vapour fraction (described in section 4.2), and the recombined hydrocarbon composition (method outlined in Aplin et al., 1999).
Fluorescence properties (i.e., λ max and red/green quotient) were measured on a Zeiss Axio Imager A2M microscope equipped with an ultraviolet (UV) light source and the Diskus-Fossil system. 20 samples from various depths were measured using Hilgers fluorescence spectroscopy, measuring the full spectrum in the visible light range (400-700 nm); this method uses two filter cubes to measure intensity of fluorescence at 600 nm (red cube) and 520-570 nm (green cube).

Raman Spectroscopy
Gas and liquid fluid inclusions were analyzed using a Renishaw inVia Raman microscope equipped with a 1,200 line/mm grating, Leica ×100 objective (N.A. = 0.85), and a motorized stage. All spectra were collected in high confocality mode to restrict Raman analysis to the smallest possible volumes with a minimum spot size of 0.76 µm and focal length of 2 µm. Acquisition times varied between 5 and 60 s per accumulation and five accumulations were collected for each measurement. The system was aligned and calibrated daily using an internal silicon reference material. Excitation via a 532 nm green laser provided~20 mW of power. Depth series analysis across fluid inclusions with a 1 µm step size (z-direction) was used to obtain the strongest Raman signal for fluid inclusion spectra. Renishaw Intelligent Background removal tool employing high order polynomial fitting was used to remove effects of fluorescence from calcite host and/or hydrocarbons within the fluid inclusions. Raman spectra were collected at relative wavenumbers of between 100 and 3,100 cm −1 . All measurements were performed at room temperature (23°C).
The Raman peak position of methane varies as a function of pressure (Fabre and Couty, 1986;Hansen et al., 2001;Lin et al., 2007;Lu et al., 2007;Zhang et al., 2015). Therefore the pressure of methane inclusions can be calculated from the shift in the C-H stretching band, v 1 , by measuring the difference between v 1 peak position in near ambient methane inclusions vs. high pressure inclusions using established relationships outlined in Lu et al. (2007) and Zhang et al. (2016). We created synthetic fluid inclusions at atmospheric pressure using fused silica capillary tubing similar to the method described by Chou et al. (2008). Briefly, a short length of deactivated fused silica tubing (0.32 mm OD) was cut and then the outer polymer coating removed before pure methane was loaded into the capillary, which was then capped using a small pieces of septa. The synthetic methane inclusions were used to calibrate the relationship between vapor pressure and methane C-H stretching Raman band position as detailed in Lu et al. (2007). A suite of three synthetic methane inclusions were analyzed giving an ambient pressure methane v 1 peak position of 2,917 cm −1 .

Stable Isotope Analysis of Carbonates
Samples for carbonate isotope analysis were milled from calcite veins and host rock using a Sherline model 4,000 lathe equipped using either a 1 mm or 0.5 mm dental drill bit. All stable isotope analyses were conducted by the Isotope Science Lab at the University of Calgary. Carbonate isotope analysis was performed via CF-IRMS on a Thermo Finnigan GasBench coupled to a Delta V isotope ratio mass spectrometer. Samples are weighed into Labco Exetainers capped and flushed with helium and then heated to 25°C at which point they were reacted with anhydrous phosphoric acid for >5hrs. CO 2 in the headspace is then sampled and analyzed six times per sample. Samples are corrected for instrumental effects using a suite of internal lab standards and normalized to VPDB using NBS-19. Values are reported in standard delta (δ) notation in ‰ VPDB, with oxygen isotope values converted to VSMOW scale (for water modeling) using the equation (Brand et al., 2014): Dolomite samples (host rock) were δ 18 O-corrected using the dolomite fractionation factor calculated after Rosenbaum and Sheppard (1986).

Fractures
Four types of fractures were observed in the core: 1) vertical to sub-vertical calcite-filled fractures with blocky calcite textures Frontiers in Earth Science | www.frontiersin.org June 2022 | Volume 10 | Article 918898 ( Figure 3A, Figure 4A); 2) highly fractured brecciated zones with a combination of vertical and horizontal fractures with a mostly blocky calcite texture ( Figure 3B); 3) horizontal bedding plane parallel fractures (beefs) composed of calcite with complex fibrous and blocky textures ( Figure 3C); and 4) sedimentfilled fractures rich in organic matter ( Figure 3D). Vertical to sub-vertical fractures are generally ≤10 mm wide, greater than 10 cm in length, but exhibit minimal vertical displacement (<10 mm) and are likely tensile fractures related to hydraulic fracturing. Bedding-plane parallel fractures are ≤10 mm thick and thought to extend over significant horizontal distances (D. Laycock, pers. comm.). Horizontal and brecciated fractures display more complex textures with evidence of multiple crack-seal episodes similar to that reported by Gasparrini et al. (2021). Vertical fractures have secondary fluid inclusions in healed fractures (see below), which when combined with complex textures in the horizontal and brecciated fractures indicates fractures have experienced multiple fluid flow and precipitation events.

Fluid Inclusions
Fluid inclusions are abundant within the calcite-filled fractures of all four samples and consist of a complex mixture of aqueous, petroleum, and gas-rich inclusions. Aqueous fluid inclusions were observed in all samples, although were rare (n = 3) in the brecciated sample (2055 m) despite the presence of abundant petroleum inclusions. Petroleum fluid inclusions were also observed in all samples, however, were rare (n = 3) in the deepest sample (2,165 m). Primary two-phase, liquid-vapour (LV), aqueous fluid inclusions typically occur as groups of large negative-crystal shaped inclusions oriented along the axis of crystal growth ( Figure 4B). These inclusions are generally isolated, do not occur in assemblages along fracture planes, and do not occur within the same fluid inclusion assemblages as petroleum inclusions ( Figures 4C,D). In addition, the absence of mixed inclusions (aqueous inclusions that contain an immiscible petroleum phase) indicates that aqueous and petroleum inclusions are not coeval. Aqueous fluid Microthermometric measurements were made on aqueous inclusions and petroleum-bearing inclusion, which includes homogenization temperatures (Th) from both aqueous and petroleum inclusions and ice melting temperatures (Tm) from aqueous inclusions as petroleum inclusions did not freeze at temperatures as low as -100°C. Eutectic temperatures (Te) were measured on several aqueous fluid inclusions to assess fluid electrolyte composition. All microthermometry results are presented in Table 1 and Figure 6. Homogenization temperatures for aqueous fluid inclusions ranged from 93°C to 149°C with a mean value of 124 ± 12°C (n = 70), and a median value of 124°C. The Th of petroleum-bearing fluid inclusions ranged from 36°C to 87°C with a mean value of 48 ± 9°C (n = 67), median value of 46°C. However, 80% of the petroleum-bearing inclusions have Th values between 40 and 50°C. There is no significant difference in Th values as a function of depth for either aqueous or petroleum-bearing fluid inclusions. The observed homogenization temperatures for both aqueous and petroleum fluid inclusions are slightly higher than the mean values reported by Gasparrini et al. (2021) for another Montney well (aqueous = 100°C; petroleum = 31°C), but both datasets have a Ice-melting temperatures (Tm) for aqueous fluid inclusions ranged from -20.0°C to -13.5°C with a mean value of -16.5 ± 2.1°C (n = 11) and exhibit minimal variability with depth ( Figure 6). These ice-melting temperatures correspond to fluid salinities of between 17.3 and 22.4 wt% equivalent NaCl, with an average salinity of 19.6 wt% NaCl equivalent. These values are nearly identical to those reported by Gasparrini et al. (2021), suggesting similar fluids were responsible for the calcite cementation at both locations. Eutectic temperatures for aqueous inclusions ranged from -53°C to -50°C (n = 7). These low eutectic temperatures are consistent with the presence of Ca 2+ or other divalent cations along with NaCl and suggests a NaCl-CaCl 2 brine composition (Bodnar, 2003).
Measurements of the vapor fraction in fluid inclusions were made using the Zeiss Zen ® Software area calculations. These are 2D measurements that are unable to account for depth and therefore have a higher degree of uncertainty. Measurement of a suite of 54 fluid inclusions indicates the average vapor fraction in aqueous inclusions is 0.045 (σ = 0.009) compared to petroleum inclusions that averaged 0.21 (σ = 0.073) vapor. Isochore-based calculations for the average aqueous inclusion molar volume (19.23 cm 3 mol −1 ) yields a fraction of vapor equal to 0.051 at 25°C, lending confidence to the accuracy of 2D vapor fraction measurements. The average petroleum inclusion vapor fraction measurement (0.21) is used in isochore modeling.

Raman Spectroscopy
All Raman spectra from fluid inclusions had significant fluorescence from the calcite host and from hydrocarbons within inclusions. Fluorescence effects were removed from Raman spectra by subtracting a polynomial curve fit to the baseline. The Raman spectra of the different fluid inclusion types are shown in Figure 7. Two phase aqueous LV fluid inclusions ( Figure 7A) have strong calcite bands from the host and no evidence of other gases such as methane or other hydrocarbons. This suggests hydrocarbons were either not present during the formation of aqueous inclusions, or immiscible; in either case the lack of CH 4 in the aqueous inclusions indicates the H 2 O+ NaCl (or the H 2 O+ NaCl + CaCl 2 ) system can be used to model aqueous inclusion trapping conditions (e.g., . Single-phase petroleum gas inclusions ( Figure 7B) are non-fluorescing and exhibit strong methane C-H stretching bands (v 1 ) at 2,905 cm −1 , smaller peaks at 2,954 and 2,890 cm −1 surround the C-H stretching peak and are interpreted as ethane or other light hydrocarbons (cf. Frezzotti et al., 2012), and strong calcite signals from the host mineral. Similar to gasrich inclusions, vapor in two-phase petroleum inclusions ( Figure 7C) display strong Raman signals in the light-hydrocarbon gas region (2,300-3,000 cm −1 ) as well as calcite bands (712 and 1,085 cm −1 ) and two broad bands from organic matter in the disordered (D) and graphite (G) band regions located at c. 1,350 cm −1 and c. 1,580 cm −1 , respectively (c.f. Henry et al., 2019). Two-phase petroleum liquids ( Figure 7D) have strong Raman signals in the organic matter D-and   (Figure 8). The carbon isotope composition of fracture-filling calcite ranges from -7.3‰ to +3.8‰ (mean of -2.6‰) and the oxygen isotope composition varies from -10.2‰ to -6.1‰ (mean of -9.0‰). One sample of a shell fragment has a δ 13 C value of -5.1‰ and δ 18 O value of -0.4‰, which is within the bounds of equilibrium with Early Triassic seawater (Veizer et al., 1999). The large difference in δ 18 O values between the shell fragment and fracture filling calcite combined with an agreement with Triassic seawater indicates there was not wholesale re-equilibrium of carbonates. Also, there is an apparent lack of exchange between host rock   Figure 9) ranging from -9.9‰ to -6.6‰. The δ 13 C values of calcite from vertical fractures vary from -7.3‰ in the Upper Montney/Doig Formation to -1.2‰ in the Lower Montney. This +6.1‰ shift in δ 13 C values with depth is associated with a +5.2‰ shift in mud-gas methane values from -49.2‰ to -44.0‰ (Figure 9). The δ 18 O values of calcite from bedding plane parallel fractures have a narrow range in values (blue in Figure 8). The δ 18 O H2O values of pore water are also plotted, alongside fracture calcite δ 18 O carb values (Figure 9).

Fluid Inclusions Isochore Modeling and Timing of Fracture Formation
Characterization of aqueous and petroleum fluid inclusions helps elucidate the physical and chemical conditions during trapping and provides a record of hydrocarbon maturation and migration within the Montney Formation. The 56°C range in homogenization temperatures from aqueous inclusions may indicate multiple generations of fluid flow, however the discrimination of these events is not readily discernable from the studied samples, and therefore they are categorized as a single generation. Petroleum inclusions occur as both primary inclusions (e.g., in the brecciated sample, 2055 m, Figures  5E-H) and as secondary assemblages within the vertical calcite-filled fractures that host the primary aqueous fluid inclusions ( Figures 5A-D). Despite occurring as primary and secondary assemblages, petroleum inclusions exhibit a narrow range in physio-chemical characteristics (e.g., homogenization temperatures, vapor fractions, fluorescence color, Raman band shifts) implying similar pressure-temperature conditions and fluid composition during entrapment. Limited range in fluorescence color may indicate either a single generation of hydrocarbons (Munz, 2001) or post-depositional alteration via phase separation or gas/water washing (Bourdet et al., 2014), however there is no evidence of gas/water washing. Low variance in the vapour-liquid ratio of petroleum fluid inclusions indicates homogeneous entrapment (Rongxi et al., 2011). This evidence strongly suggests a single generation of petroleum inclusions within these samples.
Aqueous and petroleum inclusions were not observed in the same fluid inclusion assemblages. Raman spectroscopy analysis reveals hydrocarbons are absent from aqueous inclusions, which would be expected as the presence of methane or other water-FIGURE 8 | Crossplot of δ 13 C and δ 18 O from stable isotope analysis of different carbonate phases. Fracture calcite data from this study and Gasparrini et al., (2021 FIGURE 10 | Pressure-temperature modeling for aqueous fluid inclusions. Aqueous bubblepoint curve calculated based on a measured salinity of 19.6% NaCl equiv. using the Aq_NaCl program (Bakker, 2018). Average aqueous isochore shown as dark blue line with minimum and maximum measured Th + Tm values defining light blue envelope, all values calculated using Aq_NaCl program. Green and orange lines are lithostatic and hydrostatic gradients respectively, derived from a burial history model for this well (Figure 11). Grey lithostatic and hydrostatic gradients calculated using a depth-pressure gradient of 10.18 kPa/m and a local geothermal gradient of 35°C/km (Weides and Majorowicz, 2014) based on modern conditions. Maximum burial pressure from two different burial history models are shown as grey dashed in, see text for details.
Frontiers in Earth Science | www.frontiersin.org June 2022 | Volume 10 | Article 918898 soluble organic compounds (e.g., benzene and toluene) typically occur due to interactions between water and oil in the basin (cf. Ruble et al., 1998;Gao et al., 2017). Finally, the intersection of the modeled petroleum and aqueous inclusion isochores is in excess of the pressure-temperature conditions expected for this system. Therefore, we conclude that aqueous and petroleum inclusions are not co-genetic and require separate PVT modelling.

Aqueous System
Aqueous fluid inclusion isochore and bubble point data is plotted alongside hydrostatic and lithostatic gradients based on modern geothermal gradients (35°C/km, Weides and Majorowicz, 2014) and a depth-pressure gradient of 10.18 kPa/m ( Figure 10). Modern measured reservoir pressure is between the hydrostatic and lithostatic gradient, indicating overpressure conditions currently exist within the reservoir. Also plotted are burial-history based lithostatic and hydrostatic conditions for the well, which generally agree with those calculated from the modern conditions described above. The aqueous isochore intersects with the modern formation pressure at 140°C, however this is below the hydrostatic gradient and would indicate underpressure conditions, which is unlikely. The intersection of the aqueous isochore with the hydrostatic gradient indicates a minimum trapping temperature of~142°C (i.e., ≥18°C pressure correction). Burial history for this well was derived from depositional age models and vitrinite reflectance measurements ( Figure 11) and indicates that maximum burial depths and temperatures occurred during the Laramide orogeny at~50 Ma (Watt et al., 2022;this study). Modeled trapping temperatures of~142°C indicate fracture-filling calcite formation during peak burial, as formation temperatures greater than 140°C only developed during this period. Therefore, we suggest that the vertical fractures must have formed during peak burial, associated with the Laramide orogeny, and are therefore a latediagenetic feature. Burial history models for other Montney wells, including those of Ness (2001) and Ducross et al. (2017), show deeper burial during the Laramide orogeny, potentially up to 5 km. This extends the timeframe during which formation temperatures >140°C were present, however all models indicate that the Laramide is the only event resulting in burial in excess of the 2 km that is required to produce the trapping temperatures measured in aqueous fluid inclusions. Timing of entrapment reported by Gasparrini et al. (2021) occurs prior to peak Laramide burial because their measured Th values and pressure corrections to aqueous inclusions were lower than peak burial temperatures predicted by the burial model. However, owing to the higher trapping temperatures measured in this core, initial fracture formation and calcite precipitation is suggested to occur during peak burial at approximately 50 Ma.
Eutectic and last ice melting temperatures indicate fluids with a 19 wt% NaCl-CaCl 2 brine composition, similar to those measured in aqueous fluid inclusions from Devonian and Mississippian calcite and dolomite cements (Adam, 2000;Al-Aasm and Clarke, 2004;Rivas, 2004;Davies and Smith, 2006;Al-Aasm et al., 2019). Devonian formation waters are generally composed of high-salinity calcium chloride brines, which may be the result of brine interactions with Precambrian basement rocks (Spencer, 1987). This may imply previous connectivity between these systems, however, the fluids must have been in place prior to the development of the low-permeability conditions that were present during the formation of these fractures, discussed further below.

Petroleum System
Pressure-temperature modelling of petroleum inclusions is plotted alongside the calculated phase envelope for petroleum fluids ( Figure 12). It should be noted that the phase envelope and petroleum isochore was calculated using recombined petroleum fluid compositions derived from modern produced fluids (data from operator). Therefore, it is possible that the composition of modern produced fluids is different from that within the studied petroleum fluid inclusions. However, due to the low permeability of the Montney Formation (Ghanizadeh et al., 2015), current state of overpressure conditions, and the lack direct measurement of fluid inclusion composition, we assume that the modern recombined fluid composition is reflective of the fluids within the petroleum inclusions. In addition, the PVTsim (Calsep Co.) software has been shown to produce reliable models provided the petroleum composition is at least genetically related (i.e., from the same region) to the inclusion composition (Aplin et al., 1999). The most striking feature of this modeling is the proximity of the petroleum isochore to the lithostatic gradient. Entrapment of petroleum inclusions near the lithostatic gradient indicates fluids sourced from a high-pressure system. We infer that the development of overpressure caused by hydrocarbon degradation and thermal cracking during peak burial resulted in the formation of brecciated horizons and the reactivation of previous (vertical) fracture systems. As evidenced by primary petroleum inclusions within the brecciated horizon exsolving from the host rock ( Figures 5E-H) and secondary petroleum inclusions trapped in healed fractures within the vertical fractures ( Figures 5A-D).
Furthermore, the formation of bedding-plane parallel fractures (beefs) has been shown to be associated with overpressure conditions (Rodrigues et al., 2009;Cobbold et al., 2013;Zanella et al., 2021). Beefs in shale-siltstones beds form via the development of overpressure, when pore fluid pressure exceeds lithostatic pressure at a specific depth (Cobbold et al., 2013), related to mechanical compaction (e.g., crustal shortening) and/or hydrocarbon generation/degradation (Zanella et al., 2021). The abundance of beefs in this core supports the interpretation that overpressure developed near or in excess of the lithostatic pressure and implies near synchronous formation of beef and brecciated horizons.

Raman Spectroscopy
Similarities in the Raman methane v 1 peak position between all petroleum inclusions ( Figures 7B-D) indicates these inclusions were trapped under similar pressure-temperature conditions. Indeed, the similarities in vapor fractions in petroleum fluid inclusions, the monochromatic fluorescence of these inclusions, and low variance in homogenization temperatures all suggest near synchronous formation of all petroleum inclusions independent of any depth related changes in thermal maturity.
Raman spectroscopy of petroleum inclusions reveals bathochromic downshifts to the methane C-H stretching band, v 1 , in both single-phase gas and two-phase petroleum inclusions, which indicates high pressure or high methane densities (Lin et al., 2007;Burruss et al., 2012). According to Hansen et al. (2001) who measured shifts in the methane v 1 Raman band in methane-ethane mixtures at different pressures, there is a weakening of the C-H bond, which is intensified in methane-ethane mixtures with increasing ethane contents (and would be presumably also enhanced by the presence of other hydrocarbons). This weakening of the bond caused in part by enhanced van der Waals interactions leads to charge redistribution resulting in a stronger pressure dependence on the methane v 1 band position peak (Hansen et al., 2001). The observed methane band position in this study at 2,905 cm −1 is consistent between all measured gas and petroleum inclusions and thus independent of the relative proportions of liquid and vapor (e.g., single-phase gas and two-phase) in the inclusion. This suggests that either vapor composition is the same in all fluid inclusions, or that this enhanced charge redistribution has a limit, similar to the lower limits of the pressure-methane v 1 band position relationship shown for pure methane (Fabre and Couty, 1986;Lin et al., 2005;Lu et al., 2007;Zhang et al., 2016). Notably, the methane v 1 band position of 2,905 cm −1 is consistent with methane band positions for gas hydrates (large cavity CH 4 ), and indeed we also observe the smaller 2,914 cm −1 band (small cavity CH 4 ) typically observed in methane hydrates (e.g., Sum et al., 1997). The origin of this gas hydrate-like signature is unknown.

Stable Isotopes and the Composition and Origin of Fracture Calcite
Stable isotope analysis of fracture calcite, host-rock dolomite and mud gases reveals information on the source of the fracture-filling fluids and the diagenetic processes at play. δ 13 C values of fracture calcite (vertical and horizontal) increase down core, however host rock dolomite does not (Figure 9). Dolomite formation was an early-stage diagenetic phase (Vaisblat et al., 2017;Liseroudi et al., 2020;Gasparrini et al., 2021) that formed prior to maximum burial and the onset of fracture formation. Therefore, carbon sources for the dolomitizing fluids may have been distinct from those for fracture calcite, however we only have two δ 13 C values for host rock and more work is needed to fully address this. δ 13 C values of vertical fracture calcite increases with depth concomitant, and with a similar magnitude to δ 13 C values of CH 4 in mud gas (δ 13 C values of C2 and C3 also display similar increases down hole, data not shown for clarity). Horizontalfracture calcite δ 13 C values also increase with depth, but with a greater deviation, indicating different carbon sources for cements in these two fracture types. The observed isotopic depth profile for hydrocarbon gases is representative for gradually increasing maturation of thermogenic gases (James, 1983;Schoell, 1983;Chung et al., 1988;Tilley and Muehlenbachs, 2006). Ethanepropane δ 13 C relationships are consistent with primary cracking of type II kerogen, as measured in mud gases within the WCSB (Tilley and Muehlenbach, 2006). Tandem increases in the δ 13 C values of fracture calcite suggests depth-related processes also affected fracture calcite δ 13 C composition. δ 13 C values in both vertical-and horizontal-fracture calcite increase with depth, however the magnitude of increase is greater in horizontal fracture calcite owing to the timing of calcite precipitation and associated physio-chemical conditions. From fluid-inclusion analysis, as discussed above, we know that vertical and horizontal fractures developed at different times, however, both at elevated temperatures. The negative δ 13 C values of vertical fracture calcite, in combination with negative δ 18 O values, is consistent with thermocatalytic decarboxylation and precipitation of carbonates at elevated temperatures in a diagenetic environment (Irwin et al., 1977;Klein et al., 1999;Dale et al., 2014). Increases in the δ 13 C values of mud-gas CH 4 caused by depth-related thermal maturity trends supports the depth-temperature interpretation of fracture calcite δ 13 C values.
Enhanced down-core enrichments in δ 13 C values of horizontal fracture calcite (compared with vertical fracture calcite) is likely produced by methanogenesis related to thermochemical production of CH 4 associated with hydrocarbon degradation (e.g., Hudson, 1977). CO 2 within natural gas systems that is in equilibrium with methane becomes enriched in 13 C, which is then incorporated into the carbonate pool (Hudson, 1977), and ultimately into horizontal fracture calcite reflected by increasing δ 13 C values. Hydrocarbon cracking, gas formation, and the development of overpressure at peak burial combined to enhance the incorporation of methanogenesis-equilibrated carbon in horizontal calcite producing down core enrichments in δ 13 C values compared with vertical fracture calcite.
δ 18 O values of vertical and horizontal fracture calcite are similar to host-rock dolomite values, indicating that the oxygen in fracture calcite was likely sourced from host-rock pore fluids (Figure 9). At the boundary with the Doig Formation (<1950 m), δ 18 O values of host rock and fracture calcite are elevated compared with Montney samples (Figure 9). Below this depth, δ 18 O values of both vertical-and horizontalfracture calcite do not change significantly with depth (except one sample) indicating a similar source of oxygen for all fracturefilling carbonate. The δ 18 O of pore water (data from operator) is plotted alongside fracture calcite δ 18 O and displays an increase in δ 18 O values at the Doig boundary, concomitant with δ 18 O values in fracture calcite. Further evidence that locally sourced (not migrated) fluids are responsible for precipitation of calcite in fractures. Gasparrini et al. (2021) also noted the similarities in geochemical and petrographic (e.g., cathodoluminescence) characteristics of the fracture calcite and suggested locally sourced parent fluids. To further investigate the source and relationship between fracture calcite and pore water we have modeled the oxygen isotope composition of the parent fluids to compare with the measured pore water composition.

Modelling the δ 18 O Composition of Fracture Calcite Parent Fluids
Combining aqueous inclusion trapping temperatures derived from the measured homogenization temperatures (Th) with measured calcite δ 18 O carb values, the δ 18 O H2O value of the parent fluid can be calculated (e.g., Blyth et al., 2000;Morad et al., 2010;Ardakani et al., 2013). Homogenization temperatures range from 110 to 149°C, with an average value of 123.9 ± 12.1°C; and δ 18 O carb values of vein calcite range from -9.9‰ to -6.6‰ VPDB, with an average value of -8.3 ± 1.5‰ VPDB. Using the high-temperature equilibrium fractionation equation for oxygen isotopes in carbonates (Kim and O'Neil, 1997), the oxygen isotope composition of the parent fluid is calculated to range from +5.7‰ to +12.5‰ VSMOW, with an average value of 8.6‰ VSMOW ( Figure 13). This range in modeled δ 18 O H2O values derived from aqueous inclusions agrees well with measured δ 18 O H2O values from pore waters ( Figure 13). Therefore compositionally, water within aqueous fluid inclusions appear to be in equilibrium with formation pore waters, but have clearly evolved significantly from original Triassic seawater (Veizer et al., 1999). This suggests that carbonate for calcite precipitation was likely derived locally, from pore water, and not from migrated fluids, indicating closed system behavior. This is in agreement with the evidence from the host-rock dolomite-fracture calcite δ 18 O relationship described above including previous interpretations (Gasparrini et al., 2021) and with carbon isotope records of hydrocarbon gases within the Montney (Cesar et al., 2020).
The 18 O-enriched parent fluids modeled from aqueous fluid inclusions and those measured for pore waters are similar to δ 18 O values measured in other Montney formation waters (Osselin et al., 2018;Owen et al., 2020). The NaCl-CaCl 2 brine composition (19% NaCl equiv.) and elevated δ 18 O values (compared with seawater) are similar to those observed in Devonian formation waters of the WCSB (Figure 14; Simpson, 1999). Ca-Cl brines are thought to reflect interactions of brines with Precambrian basement rocks, which provided important sources of Ca 2+ and Mg 2+ and which were potentially involved in the precipitation of dolomite cements (Spencer, 1987). However, vertical migration of fluids seems unlikely due to pervasive bedding-plane parallel permeability barriers and a lack of evidence for vertical homogenization of hydrocarbons (Watt et al., 2022). Clearly, Montney Formation pore fluids evolved from a Triassic seawater composition into Devonianlike brines prior to the development the low permeability conditions that currently exist. However, as vertical migration of fluids is unlikely we speculate that pore water brines formed via evaporative processes within the WCSB similar to those responsible for the evolution of Devonian pore waters (Connolly et al., 1990).
Fluid inclusion analysis and stable isotope geochemistry of fracture calcite and dolomite cements provide evidence that low permeability conditions existed during peak burial associated with the Laramide orogeny. The timing of when low permeability conditions developed is beyond the scope of this study, but likely occurred very early following deposition (Watt et al., 2022). Perhaps future studies on the evolution of dolomite cementation in the Montney will shed light on this intriguing question.

Implications of Closed System Behaviour
Understanding whether hydrocarbons are dominantly primary and the reservoir is self-sourced (e.g. Ardakani  Becerra et al., 2020;Euzen et al., 2021), or if hydrocarbons migrated into this system prior to the development of low permeability conditions (e.g., Sanei et al., 2015;Wood et al., 2018;Watt et al., 2022) is of high importance for reservoir evaluation. Similarities in the oxygen isotope analysis of pore waters and the modeled oxygen isotope value of calcite parent fluids implicate closed-system like behavior during the cementation of vertical fractures. Similar findings by a previous fluid inclusion study (Gasparrini et al., 2021) corroborate this interpretation, which may have significant implications for the emplacement of hydrocarbons within the reservoir. Bedparallel permeability barriers potentially implicate eastward migration of hydrocarbons from a downdip source (Laycock et al., 2021;Watt et al., 2022). However, hydrocarbon migration must predate burial during the Laramide orogeny as these hydrocarbons would need to be in place prior to thermal degradation associated with this burial event.
Alternatively, the Montney Formation could represent a hybrid reservoir with a combination of self-sourced and migrated organic matter . Nonetheless, the closed system behavior and subsequent development of overpressure associated with thermal degradation of hydrocarbons resulted in the formation of bedding-plane parallel fractures and the brecciated horizons, potentially in concert with local tectonic activity. Therefore formation of these fractures is directly related to the generation of light hydrocarbon-rich (condensate) fluids found in the petroleum inclusions (shown in this study) and currently produced from this well. In addition, the persistence of overpressure is indicative that the fracture system has done little to promote substantial fluid flow or enhance permeability in the reservoir.
Interest in the characteristics of fluid migration within fine grain siltstone reservoirs is pertinent as these types of formations are increasingly exploited for hydrocarbon production. Fluid inclusion studies of fractures within these systems play an important part in evaluating the potential role fractures play in fluid migration and enhancing permeability in these low permeability systems. The combined use of fluid inclusion microthermometry, stable isotope geochemistry, and Raman spectroscopy will be essential in the accurate evaluation of fluid flow within low permeability reservoirs and ultimately understanding their complex diagenetic history.

CONCLUSION
In this study we used fluid inclusion and geochemical techniques to better constrain the origin, timing and evolution of fluids in a complex diagenetic system. The Montney Formation has undergone multiple periods of burial and uplift related to regional tectonic activity. Associated with tectonism was the development of calcitefilled fractures resulting in an archive of fluid flow events. Differences between the characteristics of co-occurring, but not cogenetic, aqueous and petroleum fluid inclusions provide a record of these events and help elucidate the timing of FIGURE 14 | Stable isotope values from Western Canadian Sedimentary Basin (WCSB) formation waters and modeled water from Montney vertical fracture calcite. Isotope analysis from this well (yellow) and another Montney well (grey) are from pore waters extracted from cores. Montney waters from Owen et al. (2020) were derived from flowback waters. WCSB Fm. waters (Connolly et al., 1990) include Devonian, Mississippian, Jurassic, and Cretaceous formation waters, but dominantly Mesozoic formation waters. hydrocarbon maturation and migration within the basin. The main conclusions of this study are: • The formation of vertical, bedding-plane parallel and brecciated fractures occurred over a short interval during and immediately proceeding the timing of maximum burial. The lack of hydrocarbons in aqueous fluid inclusions indicates hydrocarbons were absent or in low abundance within the earliest fluids, which formed the vertical fracture system. Bedding-plane parallel and brecciated fractures postdate vertical fractures and secondary petroleum inclusions within vertical fractions represent a reactivation of this existing fracture network. • Aqueous fluid inclusions formed at elevated temperatures and pressure likely close to the period of maximum burial associated with the Laramide orogeny based on modeled trapping temperatures~142°C and burial history models. • Pressure-temperature modeling of petroleum inclusions indicates formation pressures approaching or exceeding lithostatic pressure, which is consistent with the development overpressure and the formation of bedding-plane parallel fractures. Petroleum inclusions record lower and more consistent homogenization temperatures (~45°C), have a small range in vapor fractions, and monochromatic fluorescence color under UV illumination. These similarities suggest synchronous and homogeneous formation despite occurring over a depth range of 216 m.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.