Anomalous subsurface fluid pressures, interpreted as isolated highs or lows in pore fluid potential or hydraulic head, are generally thought to signal hydraulically isolated regions. By extension, they delineate large volumes of geological media with unusually low permeability (Neuzil 1995; Ingebritsen et al. 2006) and reveal formation-scale properties not otherwise measurable. Such features have been documented in Paleozoic shales and carbonates at the Bruce nuclear site in southern Ontario, Canada, where hydraulic heads as much as 200 m below sea level and extremely large hydraulic gradients (4 m m⁻¹) between underpressured and overpressured zones have been observed (Intera Engineering Ltd 2011; Beauheim et al. 2017, 2013). The underpressured section at this particular site (henceforth the ‘Bruce site’) has been proposed to host a deep geological repository (DGR) for low- and intermediate-level nuclear waste, making it especially desirable to delineate the cause and history of the underpressure and to understand permeability there.

Analyses based on the premise of a water-saturated system at the Bruce site suggest that deglaciation-induced pore dilation can explain the underpressure and also provide evidence that the site-scale permeability in much of the section is exceptionally low (Khader & Novakowski 2014; Neuzil & Provost 2014). Likewise, extensive geochemical investigations of the site (Al et al. 2015; Clark et al. 2015; Petts et al. 2017) have indicated that the groundwater system is very old and stable, as evidenced by profiles of dissolved gases such as helium that show little migration since the Silurian Period. However, a small amount of separate phase methane gas may exist in a few horizons within the shale layers (Intera Engineering Ltd 2011; Sterling et al. 2011), so the system may not be best described as entirely water-saturated. If significant separate phase gas exists, flow and transport processes are inherently more complex and therefore introduce additional uncertainty into couplings underpinning the hydrological and hydromechanical models in these studies. Based on the salinity of the Paleozoic porewater and estimates of the thermal history of the site, any gas phase that exists there is thought to have formed a geologically long time ago (Nuclear Waste Management Organization 2011). An additional context that suggests further examination of multiphase behaviour at the Bruce site is that underpressured basins worldwide tend to contain significant gas-phase hydrocarbons (Law & Spencer 1998), perhaps pointing to a genetic link between the two.

Research in petroleum geology has produced varied conceptualizations of the relationship between gas phase and underpressures in geological basins, and some appear contradictory. For example, gas migration has been seen as leading to underpressures (e.g. Law & Dickinson 1985), while the reverse causality has also been proposed: that is, that underpressures cause gas migration (e.g. Nelson et al. 2015). The fact that such vastly different schools of thought exist may not be surprising given the complex history of many subsurface environments. Moreover, most studies of multiphase systems on large spatial scales and over geologically significant timescales have been based on geological and geochemical evidence to support qualitative descriptions of transport, rather than simulations of multiphase flow phenomena. The carefully studied and well-characterized Bruce site presents an exceptional opportunity to further investigate these issues.

In this study, we examined the role that gas phase may play at the Bruce site by using simple conceptual models of the site and its history to frame multiphase flow simulations performed with the code iTOUGH2-EOS7C (Oldenburg et al. 2004; Finsterle et al. 2017). Our goal was to discover whether mechanisms exist by which strong underpressures at the Bruce site might result from the presence and influence of gas-phase methane, and thus offer an alternative conceptual model for how the underpressures were generated. The simplified model domain, boundary conditions, initial conditions and source terms were intended to capture sufficient features of the Bruce site to address questions about the multiphase flow and transport processes that could have occurred there, and their effects on the pressure evolution of the system, without introducing unnecessary levels of complexity.

Because the Canadian Nuclear Waste Management Organization (NWMO) has proposed a Deep Geological Repository (DGR) at the Bruce site, extensive and careful characterization of the local and regional geology, hydrology, and geochemistry has been carried out there. Located on the eastern flank of the Michigan Basin, the proposed DGR would be at a depth of c. 680 m BGS (metres below ground surface) within the Ordovician Cobourg Formation, an argillaceous limestone overlain by more than 200 m of shale. As Figure 1 shows, straddle packer tests in vertical and inclined boreholes at the site have revealed extremely low hydraulic conductivity parallel to bedding (K_h), and a distinct underpressure within the Ordovician-aged shale and limestone sequence between the depths of c. 395 and c. 850 m BGS. The black lines in Figure 1a show formation-average values of K_h calculated by Walsh (2011), based on more extensive measurements reported by Roberts et al. (2011). Because the more permeable units above and below the low-permeability region, namely the Guelph and Cambrian sandstone aquifers, respectively, both display overpressures, it has been suggested that the underpressure is superimposed upon an otherwise mildly overpressured regime (Neuzil & Provost 2014). A coarse approximation of the nominal overpressure is indicated by the blue line in Figure 1b.

Hydraulic conductivities (Fig. 1a) in strata that are underpressured (Fig. 1b) are extremely low (the harmonic mean indicated by the red line in Figure 1a is equivalent to a permeability of c. 2 × 10⁻²¹ m²; see also the hydraulic conductivity of the Cambrian aquifer for comparison). The underpressure is pronounced, reaching c. 5 MPa below brine hydrostatic, with density-corrected environmental heads as much as c. 200 m below sea level (Intera Engineering Ltd 2011). The underpressure is remarkably consistent among the four vertical boreholes, which are spread across an area of about 1.6 km² surrounding the perimeter of the proposed DGR. Head gradients indicated by the pressure data are unusually steep, especially near the bottom of the underpressured zone, and became more pronounced during several years of monitoring (Beauheim et al. 2017). Because equipment failure in such an extremely anomalous system would almost certainly drive pressure observations back towards hydrostatic, the fact that the measurements have instead deviated further from it during the measurement period indicates that the measurements are still approaching the pre-drilling conditions. The coloured horizontal labels in Figure 1a give the names of the layers in the simplified model domain developed for this study, as described in the next section.

Permeability anisotropy presumably exists at core to formation scales, but is not well characterized. Differences in measured K_h in a single horizon, moreover, are often comparable to the differences seen between geological units. Rather than estimating vertical hydraulic conductivity, K_v, via scaling calculations that are based on assumed anisotropy ratios and thus invoke additional uncertainty, we used the originally measured K_h values to estimate the permeabilities of the layers in our model, as shown by the red and blue vertical lines and the green asterisk in Figure 1a. These values are likely to represent upper-bound estimates, but sensitivity analysis simulations (not shown herein) revealed that decreasing the permeabilities by an order of magnitude had a relatively small effect on the results of this study.

Different lines of evidence have raised the possibility that gas-phase methane is present within the underpressured Ordovician rocks at the Bruce site. First, measurements on core samples by different laboratories indicated the presence of a gas phase at saturations that were relatively consistent: for example, at an average value of about 10–15% in the Ordovician units (Sterling et al. 2011). These findings were based on two different petrophysical testing techniques for measuring liquid saturation (the standard Dean Stark extraction method and the nuclear magnetic resonance (NMR) method), on geochemical calculations indicating dissolved gas concentrations in porewater above solubility limits, and on calculations indicating that the volume of gas within the core samples exceeded the volume of microcracks that probably formed during coring and sample preparation procedures. How indicative these laboratory techniques are for the extremely low-porosity, low-permeability and low-gas saturation conditions at the Bruce site, however, is unclear. Dilation of the cores and evaporative loss of porewater are difficult to account for, and borehole indications of methane gas occurred where organic matter contents are highest, suggesting that the methane was initially adsorbed (Neuzil & Provost 2014).

Second, borehole hydraulic tests indicated higher than expected test zone compressibility (C_tz) in some horizons, including within the shale, as shown in Figure 2. These observations may be most readily explained by the presence of gas in the borehole (Beauheim et al. 2013); however, it is unclear whether the gas was present before drilling or simply exsolved from porewater in the surrounding rock during the drilling process. If gas-phase methane is naturally present in the rock, it may be limited to a few stratigraphic horizons where the strongest suggestions of gas-phase methane have been found (Intera Engineering Ltd 2011; Sterling et al. 2011). As shown by the orange lines in Figure 2, C_tz measurements statistically higher than the rest were found in three zones that occupied depth ranges of c. 270–330, c. 430–480 and c. 570–660 m BGS. The grey lines show the means of the anomalously high C_tz measurements within these stratigraphic groupings.

Various researchers have used palaeohydrogeological modelling to test different hypotheses regarding potential causes for the underpressure at the Bruce site. While Neuzil & Provost (2014) showed that glacial loading cycles could plausibly create the observed underpressure in a water-saturated system, other studies have included a gas phase in their analyses. Normani & Sykes (2012) summarized a set of hydromechanical simulations that represented the evolution of the site over a 120 kyr-long period including the Pleistocene ice-sheet advance and retreat, in which loading efficiencies and storage coefficients were adjusted to represent the presence of a gas phase. Although underpressures did result, they were much less pronounced than those observed, and multiphase flow physics were not explicitly accounted for. Previously, Sykes et al. (2011) analysed effects of multiphase flow at the site and were able to match the observed underpressure, for certain specific assumptions and conditions. In one of their simulations, high gas saturations and water pressures of zero (underpressures) were set as initial conditions across most of the Ordovician section. In the other, transient gas generation was simulated from an initially hydrostatic state, but underpressures only developed when the Millington & Quirk (1961) model – one with unknown applicability to the types of rocks encountered at the site – was used to represent multiphase diffusion. Questions thus remain about the transient evolution of the gas, the initial development of the underpressure and the connection between the two.

We suspected that gas-phase formation and concurrent migration with co-resident groundwater could not entirely explain the underpressure observed at the Bruce site; however, we considered it desirable to further explore this possibility. Our hypothesis was that the introduction of a second fluid (gas) phase into previously water-saturated rock would not act to decrease the pressure of either fluid. Instead, as long as the gas-phase saturation remains realistically small, it could reasonably be expected that the water phase would remain relatively connected and mobile, preventing significant perturbations to its pressure. These postulations pertain only to water pressure perturbations caused purely by the presence of a gas phase, and thus do not preclude formation of pressure anomalies by other means.

A simple multiphase flow model was used to consider this issue and, more specifically, to test the effects of different rock properties, initial and boundary conditions on gas migration, and pressure evolution processes over geological timescales in a low-permeability domain like that at the Bruce site. Considering the aforementioned heterogeneity, and the fact that vertical head gradients at the site far exceed horizontal, we hypothesize that flow can be adequately represented as vertical. This allowed the use of a 1D model domain, and also made it easier to track multiphase flow processes that may have occurred at the Bruce site during the simulations. The model domain represented the rock section between 230 and 860 m BGS, using a total of 63 grid blocks that were each 10 m high.

In the absence of other influences, such as topographically driven flow, pressure in the liquid (water) phase should remain near hydrostatic and the gas pressure should therefore be greater than hydrostatic by an amount approximately equal to the saturation-relevant capillary pressure. Small gas saturations are to be expected in low-permeability rocks based on the high resistance to gas penetration provided by very high capillary entry pressures in characteristically small pores. These postulations are formulated in the context of saturation-dependent capillary pressure and relative permeability relationships that form the basis for conventional multiphase flow theory, which has been developed and tested extensively, albeit mainly for relatively high-permeability granular media such as aquifers and conventional petroleum reservoirs. Use of these relationships in the Bruce site rocks, which have exceptionally low porosity and permeability, is thus a potential limitation that should be borne in mind.

The model domain described above was implemented in simulations performed with the inverse Transport Of Unsaturated Groundwater and Heat (iTOUGH2) finite-difference multiphase flow simulator (Finsterle et al. 2017), coupled with the thermodynamic fluid property estimation module EOS7C (Oldenburg et al. 2004).The iTOUGH2 code solves mass conservation equations for each chemical component of the system (e.g. methane and water), while the EOS7C module determines the partitioning of components among fluid phases (e.g. liquid and gas). The coupled modelling tool represents flow, transport and fluid interaction processes via constitutive laws, namely the multiphase extension of Darcy’s law for advection, Fick’s law for diffusion, and the Peng–Robinson equation of state for phase partitioning of chemical components and determination of fluid properties (i.e. density, viscosity and enthalpy). The code neglects mechanical dispersion, but omitting this phenomenon was assumed to be acceptable because the system is considered to be diffusion-dominated. For this analysis, the inverse capabilities of iTOUGH2 were not used and the code was implemented in forward mode, with system parameters specified a priori.

The simulations invoke a number of assumptions, both explicit and implicit. First, although other constituents are present in situ, the problem was simplified to a pure water–methane system with only gas and liquid phases, and chemical component partitioning was modelled using the local equilibrium assumption, which is justified due to the long timescales of the analysis. Second, it was assumed that the mineralogical compositions of the rocks did not affect the properties of the fluids in any way (i.e. the rocks were treated as non-reactive, non-sorbing and of uniform, non-hysteretic wettability). Furthermore, the simulations were performed under isothermal conditions because we considered thermal effects to be second order. Thus, thermally induced fluid expansion/contraction and exsolution/dissolution processes were not accounted for. Gas-evolution processes (including exsolution, dissolution and free-phase migration) were limited to those coupled with pressure and concentration conditions, allowing for elucidation of their fundamental relationship with pressure evolution.

It was also assumed that the equation developed by van Genuchten (1980) could be used to describe the capillary pressure–saturation relationship for the rocks at the Bruce site, and that the relative permeability–saturation relationships developed by Mualem (1976) and Corey (1954) were applicable. These constitutive equations are shown in equations (1–3):

Here, P_cap is the capillary pressure, P₀ is the entry pressure, k_rl is the liquid-phase relative permeability, k_rg is the gas-phase relative permeability and λ is a fitting parameter that defines the shapes of the curves. In equations (1–3), S* and $S^$ are two different normalized saturations, defined by:

where S_l is the liquid-phase (water) saturation, S_ls is the saturated liquid saturation, S_lr is the residual liquid saturation and S_gr is the residual gas saturation. As already noted, use of equations (1–5) is a significant assumption because of the exceptionally small porosity, permeability and pore throat size of rocks represented in the simulations. The code requires a maximum capillary pressure constraint but the simulated scenarios never approached this limit, so it did not influence the results.

Rock units were characterized based on multiphase flow characteristics of the strata at the Bruce site as measured in previous studies (Jackson & Wigston 2010; Whitney & Lee 2010; Raven & Jackson 2011). Mercury porosimetry characterizations of capillary pressure, compiled in a single plot, are shown in Figure 3. Inspection reveals three stratigraphic zones with distinct capillary pressure–saturation characteristics, represented by the three colours in the figure. The zone that showed the highest capillary pressures is labelled ‘tight’, while those with more moderate capillary pressures are ‘medium’ and the Cambrian aquifer has relatively low capillary pressures.

To represent multiphase flow properties of the three regions of the model, the van Genuchten capillary pressure–saturation function was fitted to each lumped dataset, as shown by the solid lines in Figure 3. Beyond P₀ and λ, which were the two parameters used to fit the van Genuchten formula to the capillary pressure data, the only additional parameter needed to define the relative permeability curves for these lumped groups was the residual gas saturation, S_gr. Values for this parameter were taken as the average of those reported by Calder (2011) for samples that came from within each stratigraphic region. These S_gr values were simply subtracted from 1 to yield values for S_ls that were used in the relative permeability functions, which effectively allowed for an immobile gas phase to develop and persist during the simulations. In all of the capillary pressure functions, S_ls was set to 1 to be consistent with the fully water-saturated initial conditions used in the simulations, and to allow the gas phase to dissolve into undersaturated water (a necessary precursor to diffusive flux through the mostly water-saturated system).

Permeabilities for each region of the model domain were set at the harmonic mean of the field-measured permeabilities in that section, as described above (see Fig. 1); and more basic physical parameters such as porosity were taken as simple averages of the laboratory-measured values. In the simulations, the diffusive flux of component κ (in this case, methane or water) in the fluid-phase β (in this case, liquid or gas) was calculated from the following adaptation of Fick’s law (Pruess et al. 2012):

where τ₀ and τ_β are the porous medium- and saturation-dependent tortuosity factors, respectively, ϕ is porosity, ρ_β is the fluid-phase density, $dβκ$ is the diffusion coefficient of the component in the bulk fluid phase and $Xβκ$ is the mass fraction of the component in the fluid phase. The diffusion coefficient for methane in water was globally set to 1.47 × 10⁻⁹ m² s⁻¹ (Maharajh & Walkley 1973), τ₀ was calculated for each region of the model domain based on effective diffusion coefficients and porosities reported by Al et al. (2015), and the code internally set τ_β equal to k_rβ(S_β), as calculated by equations (2) and (3), in each grid block at each time step.

The resulting system parameters are summarized in Table 1. The simplified model geometry, as noted earlier, is in Figure 1. The multiphase diffusion model described by equation (6) is perhaps more appropriate for the very-low-porosity and low-permeability rocks at the Bruce site than the Millington & Quirk (1961) model, which previous studies have used, because in this type of system pockets of gas phase are likely to be small and isolated. This disconnected distribution of free-phase gas, along with the relatively low solubility of methane compared to other gases such as CO₂, should theoretically lead to slower diffusion, which the relative permeability model would reflect, rather than the enhanced diffusion process that the Millington & Quirk model would predict.

Six different simulations were conducted to test the effects of various factors on gas phase and pressure evolution through geological time in the Bruce-like system. The details of the simulations, including geological configurations and associated adjustments to the rock properties, as well as initial conditions, boundary conditions and source terms, are summarized in Table 2. Model input and output files were made publically available as part of a USGS data release (Plampin & Neuzil 2018). The simulation scenarios and their motivations are briefly described below, while their results are considered together in the next section to facilitate comparison.

Simulation 1 tested the effects of transient gas-phase methane generation and migration on the evolution of pressure within the simple Bruce-like system through geological time, using the ‘homogeneous’ geological configuration determined from the analysis described above. This ‘base-case’ simulation, and all the others described below except for simulation 6, began with a fully water-saturated system and a hydrostatic water pressure gradient, then gas-phase methane was introduced uniformly across the entire tight region at a total rate of 2.3 × 10⁻¹¹ kg s⁻¹ for a period of 100 kyr. Although this gas generation scheme caused only small gas saturations to develop (c. 1%, as shown in a later section), it represents the largest perturbation that did not cause gas pressure to significantly exceed the least principal stress, as approximated by the vertical stress. This was done because of the possibility that higher gas pressures could have fractured the rock and permitted gas to escape, whereas there is no evidence from cores and borehole logs for extensive fracturing at the site (Engelder 2011; Intera Engineering Ltd 2011). This constraint was conservative, because effective stress in systems with multiple fluid phases is prorated to the fluid saturations (Bishop & Blight 1963), and the phase with the higher pressure (i.e. the one that would cause fracturing if sufficiently pressurized – the gas phase in this case) was effectively limited to very small saturations in the simulations.

The simulations were run well beyond the transient gas generation period to allow the fluids substantial time (e.g. 10 myr in simulation 1) to redistribute and interact. This general simulation scheme is qualitatively consistent with geochemical evidence from the Bruce site (Intera Engineering Ltd 2011), which suggests that methane could have been produced biologically during the early Paleozoic Era, when the Ordovician sediments were still geologically young (Nuclear Waste Management Organization 2011). Furthermore, geological evidence suggests that up to 1.5 km of overlying material could have subsequently eroded from the surface at the site. The resulting decrease in stress could have led to further gas-phase exsolution and/or expansion. This study does not address these processes specifically, but rather looks at how the presence of gas-phase methane could affect the system during its more recent geological history.

The effects of geological heterogeneity on gas evolution and pressure distribution were considered in simulations 2–4. Simulations 2 and 3 incorporated ‘mild’ heterogeneity, while simulation 4 incorporated ‘extreme’ heterogeneity. In all of these simulations, relatively well-determined rock properties (Intera Engineering Ltd 2011) including density, porosity and permeability were left unchanged from the base case (simulation 1), and only the more poorly-constrained multiphase flow parameters were adjusted within sub-regions of the domain. These parameter adjustments were motivated and informed by the variability in the laboratory-measured capillary pressure data from the core samples.

Figure 4 shows capillary pressure curves used to represent mild heterogeneity, superimposed on the original data. The van Genuchten parameter 1/P₀ was varied by a factor of 2 in either direction from the original value (increased in simulation 2 and decreased in simulation 3) in order to capture the spread among samples in the regions of the model domain. This resulted in new capillary pressure–saturation curves labelled as ‘lower bound’ (used in simulation 2) and ‘upper bound’ (used in simulation 3). The original fitted curves from Figure 3 are shown as a ‘lumped fit’, like before. The adjusted parameters were applied within regions of the model domain where borehole hydraulic tests showed relatively high C_tz (see Fig. 2), because these zones were thought to potentially contain higher gas saturations, and entry capillary pressure, as quantified in iTOUGH2 by 1/P₀, has a large influence on gas evolution in porous media (Plampin et al. 2014). Two of these regions were in the tight region of the simplified domain, while one was within the medium region. Parameter adjustments were therefore made in only these zones, and the properties of the Cambrian aquifer were not changed.

More pronounced heterogeneity was incorporated into simulation 4, as informed by data from a few ‘outlier’ samples omitted from the previously described parameter estimation process. Figure 5 shows these data as black symbols, superimposed on the earlier data and labelled in the legend with the depths from which the samples came. The van Genuchten formula was fitted to each of these sets of data, as shown by the black curves, and each set of outlier parameters was applied to a single 10 m-high grid cell at the corresponding depth, representing an isolated subunit within the tight stratigraphic section.

Various phenomena, such as cycles of glacial ice cover, uplift and erosion, could cause pore fluid pressure in the aquifers at the Bruce site to change during the period simulated. Although it was outside the scope of this study to systematically examine such effects, we did consider a single scenario with changing aquifer pressures in simulation 5. This simulation used the homogeneous domain, methane gas generation history and initial water-saturated state of simulation 1, but applied transient boundary conditions. Specifically, the simulation began with a uniform overpressure that did not exceed lithostatic, namely a hydrostatic profile with a uniform c. 2.7 MPa overpressure. After the 100 kyr-long gas introduction period, water pressure in the top grid cell was decreased to the hydrostatic value for that depth. After 100 kyr of fluid redistribution (i.e. at t = 200 kyr), pressure in the bottom grid cell was decreased to hydrostatic for that depth, and the system was then allowed to evolve until a final time of t = 300 kyr. This scenario was not intended to test any particular hypothesis about the physical history of the Bruce site, but rather to elucidate the fundamental effects that transient boundary conditions, potentially imposed on the system during hypothetical changes to its surroundings (e.g. glacial loading cycles), may have on the processes under investigation.

Because the state of the system at the start of our simulations is unknown, simulation 6 was performed with an alternative set of initial conditions. Similar to previous investigations of the site (Sykes et al. 2011; Normani & Sykes 2012), this simulation incorporated a uniform gas saturation (S_g) of 0.01 within the tight region at t = 0. To ensure that pressure would not exceed lithostatic, and to remain consistent with the theoretical result of a small gas-phase methane density, a vertical gas pressure gradient was applied throughout the tight region, with a magnitude equal to hydrostatic water pressure in the water-saturated Cambrian aquifer where capillary pressure was zero. This initial state implicitly required underpressures in the water phase throughout the tight section due to the relatively small difference between the vertical ‘gas static’ gas-pressure profile and hydrostatic, compared to the very high capillary pressures caused by even a small gas saturation in the tight media. Since a gas phase was present initially, no methane source term was incorporated into simulation 6. Instead, the system was simply allowed to evolve from its initial state for 100 kyr.

Figure 6 shows the results of the simulations performed for this study. Each row of graphs shows, from left to right, the evolution of water saturation, aqueous methane (CH₄) concentration, gas pressure and water pressure throughout the course of at least one simulation. Showing the results of all of the simulations together in this way facilitates comparison and shows the effects that various factors had on system evolution. Figure 6a shows the results of the base-case simulation (simulation 1). Figure 6b shows the effects of ‘mild’ geological heterogeneity by comparing simulation 2 (for lower-bound capillary pressure curves applied in regions of high C_tz) with simulation 3 (for upper-bound capillary pressure curves applied in regions of high C_tz). Figure 6c shows the results of simulation 4 (using curves fitted to the outlier samples applied at the depths from which those samples came) and thus elucidates the effects of ‘extreme’ geological heterogeneity by comparison with the results of simulation 1 (base case) in Figure 6a. Figure 6d shows the effects of the transient boundary conditions applied in simulation 5, and Figure 6e shows the effects of the gas-phase saturation and vertical gas-pressure profile incorporated as initial conditions in simulation 6.

Figure 6a–c reveals that the tested geological heterogeneities had relatively small and short-lived effects on water pressure, with the latter at or close to hydrostatic throughout the simulations. Furthermore, while the heterogeneities did lead to significant differences in gas-phase formation and accumulation at early times, accumulated gas-phase methane redissolved and formed very similar water saturation (and, consequently, gas pressure) distributions before 10 myr in these simulations.

Figure 6d shows that the transient boundary conditions applied in simulation 5 did not lead to water-pressure profiles that deviated greatly from steady-state linear head gradients. The gradients in the aquifers remained very close to hydrostatic, and slight overpressures were observed within the tight section. Only simulation 6, shown in Figure 6e, produced significant deviations, as underpressures, from hydrostatic. As in previous modelling studies (e.g. Sykes et al. 2011), however, underpressure resulted from an imposed initial condition, shown by the orange lines. We were unable to attain significant underpressures in the course of any simulation by introducing a gas phase without generating gas pressures well in excess of lithostatic stress. Another important consideration is that the underpressure in simulation 6 dissipated relatively quickly; that is, the water pressure returned to a near-hydrostatic profile within 100 kyr.

In all of the simulations, minimal redistribution of the gas phase occurred because the small gas saturations kept the gas-phase relative permeability small. The methane instead dissolved into the water and then slowly diffused out of the tight region into the bounding aquifers. It is not expected that this conclusion would change even if the simulations had involved higher gas saturations (i.e. closer to the values of c. 10–15% that were estimated from core testing), because free-phase gas migration would have still been strongly limited by the very low absolute permeabilities of the rocks. We conclude that the presence of a gas phase at the saturations considered here is unlikely to completely explain the underpressures at the Bruce site.

Strong underpressures and elevated levels of methane in Paleozoic sedimentary rocks at the Bruce site in southern Ontario, Canada raise the possibility that a gas phase is present in situ and has a causal link to the underpressure. A gas phase in addition to liquid water complicates the analysis of flow and transport in the subsurface environment. In addition, applying existing multiphase flow modelling tools to environments like the Bruce site, where the porosities and permeabilities are much lower than in systems typically studied, should be done recognizing its limitations.

With these considerations in mind, we analysed long-term evolution of a water–methane system in a simplified domain with attributes of the Bruce site. We sought to examine whether underpressures in this region might be generated without externally-imposed mechanical deformation and thus did not consider stress changes imposed by erosion or cycles of glaciation. Preliminary efforts to study the latter have been reported by Walsh et al. (2012). A highly simplified 1D representation of the subsurface was used here to aid interpretation of the simulations. Although ample data are available to allow refinement of the hydrostratigraphy within the sedimentary sequence (Calder 2011; Intera Engineering Ltd 2011; Roberts et al. 2011), we do not expect that incorporating such refinements would substantially change our findings. Likewise, the simulations included only pure water and methane gas, but the salinity of brine at the site (Nuclear Waste Management Organization 2011) could have significant effects on gas solubility (Duan & Mao 2006). Further analyses of the effects of brine-relevant gas solubility on multiphase methane and pressure at the site would therefore also be beneficial (Akhbari & Hesse 2017). Due to the common co-occurrence of gas phase and underpressures in other sedimentary basins around the world (Vinard 1988; Corbet & Bethke 1992; Law & Spencer 1998), examples of which can be found in both glaciated (Masters 1979) and non-glaciated (Davis 1984) regions, improved understanding of the Bruce site may inform efforts to understand those systems as well.

The simulations tested the effects of various types of geological heterogeneities, initial conditions and boundary conditions on the evolution of fluid pressures over times of the order of 10⁷ years. Such time periods are modest geologically, but much longer than typically considered in engineering-based simulations of multiphase systems. Recognizing the limitations of the analysis, the results do not suggest that multiphase flow processes alone can generate the underpressure observed at the Bruce site. The only simulation that produced noticeable underpressure began with underpressure as an initial condition, a result of the chosen initial gas saturation and gas-pressure distribution. Even in this relatively arbitrary case, the anomalous pressure dissipated within 10⁵ years. Given that the permeabilities used in the model may be upper-bound estimates, this period could feasibly last perhaps an order of magnitude longer, but the artificial initial conditions would still require further explanation.

Although laboratory analyses of core samples have suggested non-zero gas-phase saturations throughout the Ordovician section at the Bruce site (Intera Engineering Ltd 2011), it is possible that methane may not have existed as a gas phase prior to drilling disturbance, or that it may exist only in a few isolated places such that the behaviour of the overall system is that of a water-saturated system. In either case, analyses have shown that glacial load cycling is a plausible cause of the underpressures (Khader & Novakowski 2014; Neuzil & Provost 2014). Alternatively, more complex, coupled physicochemical processes may be necessary to understand the origin and persistence of underpressures in low-permeability geological systems like the Bruce site that contain a gas phase. Potential future studies may investigate the influence of alternative representations of multiphase diffusion, salinity, sorption and capillary hysteresis on gas-phase evolution, as well as these factors in tandem with coupled hydromechanical processes in which a gas phase is explicitly represented.

The extremely low porosities and permeabilities of the Ordovician sediments underlying the Bruce site near Tiverton, Ontario, Canada, as well as the lack of pervasive fractures therein, have led to a very stable groundwater flow system in which diffusion-dominated solute transport processes occur at extremely slow rates over geological timescales. The findings of the present work are consistent with these important characteristics of the site.

The authors thank Mark Jensen for his support and critique of this work, and Dylan Luhowy for providing the most recent pressure data. Robert Walsh, Alden Provost and one anonymous reviewer are acknowledged for their insightful critiques of this work. Data used herein were not gathered by the United States Geological Survey. Any use of trade, product or firm names in this publication is for descriptive purposes only and does not imply endorsement by the United States Government.

The Canadian Nuclear Waste Management Organization is acknowledged for providing funding to C.E. Neuzil for this work through Technical Assistance Agreement 14HWTAACE00002.