Velocity modeling of supercritical pore fluids through porous media under reservoir conditions with applications for petroleum secondary migration and carbon sequestration plumes

Secondary migration of pore fluids through permeable carrier beds or along higher permeability fault and fracture systems is traditionally modeled using time-invariant transport equations, which assume instantaneous fluid flow. This time step is approached differently for modeling fluid transport processes in petrophysics and petroleum engineering occurring at production timeframes compared to secondary migration modeling for basin-scale analysis of petroleum systems over geologic timeframes (Figure 1). As an example of this differential treatment of time steps among disciplines, governing equations for natural gas transport in shale reservoirs are described by four scales of distance and time. These include molecular-collision fluid dynamics within the kerogen, microscale flow within the low-permeability matrix, grain-scale flow within the microfracture network, and mesoscale flow in the formation (Javadpour et al., 2007; Wang and Reed, 2009; Ozkan et al., 2010; Wei et al., 2018).

Carruthers (2003) determines that several input variables required for fluid transport modeling at reservoir scales and production time horizons are irrelevant for flow modeling of secondary migration at basin scales. The appropriate fluid transport equations for secondary migration in basin-scale modeling are described by buoyancy forces, capillary pressure, statistically defined invasion percolation, and a time-invariant form of Darcy’s law (Hunt, 1996; Hantschel and Kauerauf, 2009). A harmonic average of buoyancy and capillary forces is appropriate to consider for secondary migration occurring vertically across stratigraphic layers (Schowalter, 1979); however, velocity is time invariant and assumed to be instantaneous.

In this study, a time-dependent form of Darcy’s law was used to model fluid flow velocity within a low-permeability carrier bed. Fluid behaviors at reservoir conditions were determined from equations of state. The timeframe of secondary migration from source to reservoir is important in petroleum systems modeling to determine with respect to petroleum generation, fault movement, and seal formation. In addition, quantifying pore fluid migration velocities has direct applications in evaluating formations for geologic sequestration and long-term storage of carbon dioxide greenhouse gas.

In contrast to traditional methods for modeling secondary migration, which assume time invariance of fluid transport, this study’s numerical modeling objectives were fourfold: (1) quantify secondary migration velocities for supercritical fluids through porous media under reservoir pressure and temperature conditions, (2) use these velocities to calculate secondary migration distances within a given carrier bed, (3) use secondary migration velocities to calculate the transit time of these fluids, and (4) evaluate long-term storage and containment of sequestered fluids based on permeability-based flow units.

Numerical models for methane and carbon dioxide were developed to study secondary migration velocity over a broad range of reservoir pressures (0–200 MPa [0–29,000 psi]). A reservoir temperature of 177°C ± 6°C (350°F ± 20°F) was chosen for case study 1 involving methane , and a reservoir temperature of 93°C ± 6°C (200°F ± 20°F) was chosen for case studies 2 and 3 involving CO₂ migration, based on values available in research studies (Burke, 2011; Wang et al., 2013). The spectrum of porosities (⁠$ϕ$ = 0%–95%) and 15 permeability magnitudes (⁠$κ$ = 10 D–0.1 nD) encompass a substantial range of possible rock property combinations encountered in the sedimentary sequence. Three case studies demonstrate the application and interpretation of these modeling results specifically for (1) quantifying migration velocities of methane within a low-permeability carrier bed exhibiting representative rock properties of the Upper Jurassic Haynesville Formation of the U.S. Gulf Coast petroleum system, (2) defining permeability-based flow units of saline formations that are potentially targeted for geologic carbon sequestration, and (3) estimating long-term migration distances and containment timescales of geologically sequestered carbon dioxide plumes. Given the broad range of properties represented in these models, velocity results are applicable for an array of methane-bearing petroleum systems and geologic sequestration targets.

The fluid transport methodology introduced in this study is summarized in Figure 2. In continuum fluid dynamics, the Navier-Stokes momentum equations are the constitutive relations for describing the flow of a Newtonian fluid. Fluids typically encountered in the sedimentary sequence adhere to Newtonian fluid behavior in which viscosity remains constant with the application of shear stress; a rare exception is emulsified bitumen (Bazyleva et al., 2010). Darcy’s law is a general solution of the Navier-Stokes equations that are appropriate for describing the laminar flow of a viscous Newtonian fluid through a homogeneous isotropic media under a pressure differential. A time-dependent form of Darcy’s law was used to quantify flow vectors in this study.

Constitutive laws used to describe thermophysical properties of supercritical fluids under reservoir pressure and temperature conditions are discussed in the following section. Matrix porosity and permeability are used to condition the momentum equations. Model outputs quantify migration timescales. By using a 1.0 km distance, interstitial pore fluid velocities can be obtained.

Constitutive laws describing thermophysical properties as a function of pressure and temperature are used to condition the Navier Stokes and the continuity momentum equations (Figure 2). Thermophysical properties for fluid density, fluid viscosity, fluid compressional-wave (P-wave) velocity, and fluid compressibility shown for methane (Figure 3a–3d) and for CO₂ (Figure 3e–3h) are obtained using the reference fluid thermodynamic and transport properties database (NIST, 2019). Equations of state for supercritical methane (Younglove and Ely, 1987; Setzmann and Wagner, 1991; Richter and McLinden, 2014) are used to obtain thermophysical properties based on the reference fluid thermodynamic and transport properties database (NIST, 2019) for pressures ranging from 0MPa to 200 MPa (0–29,000 psi) at isothermal conditions of 350°F (177°C). Similarly, equations of state for supercritical CO₂ (Span and Wagner, 1996; Lemmon et al., 2011) are used to obtain CO₂ thermophysical properties using the thermodynamic database (NIST, 2019) for pressures ranging from 0 MPa to 200 MPa (0–29,000 psi) at isothermal conditions of 93°C (200°F). This broad pressure range is inclusive of pressure conditions likely encountered in the sedimentary sequence, considering that the deepest producing wells reach 9144 m (30,000 ft; Dyman and Cook, 1998) and the standard definition of a 1.0 psi/ft lithostatic pressure gradient at depths of 9144 m (30,000 ft) yields 206.8 MPa (30,000 psi/ft).

Fluid densities for supercritical methane and CO₂ (Figure 3a and 3e) rapidly increase with pressure. With increasing pressure, the viscosity of methane exhibits almost linear behavior (Figure 3b), whereas CO₂ viscosity approaches asymptotic behavior approximately 25 MPa (Figure 3f). The fluid P-wave velocity of supercritical methane increases linearly as a function of pressure (Figure 3c), whereas the velocity of CO₂ rapidly increases with pressure conditions (Figure 3g). Fluid compressibility for both fluids is strongly influenced by increasing pressure up to approximately 25 MPa, at which point compressibility attains a minimum value (Figure 3d and 3h).

Thermophysical properties of methane and CO₂ vary as a function of volume, pressure, and temperature. Calculations of thermophysical properties from NIST (2019) require an assumption to be made for either isochoric, isopressure, or isothermal conditions. Pore pressure increases in the rocks from, for example, buoyancy or volumetric expansion as fluid propagates, are taken into account by the pressure component in the Darcy equation. Pressure conditions are expected to have a much greater influence on the thermophysical properties of interest in this migration model (i.e., fluid density, viscosity P-wave velocity, and compressibility) than the small temperature differentials associated with geothermal gradients during migration within a carrier bed. Thus, a wide range of pressures are incorporated into the models by making the reasonable assumption of isothermal conditions. The “Sensitivity analysis” section discusses the implications of this assumption on the thermophysical properties and resultant migration velocities.

Secondary migration velocities are calculated for a substantial range of rock and fluid properties and are shown for 15 orders of magnitude in permeability for a wide range of porosities from 0% to 95% (Figure 4). Figure 4a shows that carrier beds with permeabilities in the Darcy (D) range correspond to migration velocities from 3 × 10⁻⁴ to 3 × 10⁻¹ m/s. Carrier beds with millidarcy (mD) permeabilities (Figure 4b) correspond to migration velocities from 3 × 10⁻⁷ to 3 × 10⁻⁴ m/s. Migration velocities within carrier beds measured in microdarcy (μD) permeabilities (Figure 4c) fall within the 3 × 10⁻¹⁰ to 3 × 10⁻⁷ m/s velocity range. Figure 4d shows that rock properties in the nanodarcys (nD) range correspond to migration velocities from 3 × 10⁻¹³ to 3 × 10⁻¹⁰ m/s. For completeness, permeabilities in the picodarcy range (not shown) represent the extreme lower bound in this model and correspond to migration velocities slower than 3 × 10⁻¹³ m/s, which is equivalent to a timeframe of billions of years, exceeding the tectonic cycle of most sedimentary basins. Depending on the formation permeability of interest, Figure 4 can be used to research and estimate the migration velocity by any given porosity value.

By normalizing migration through a 1.0 km carrier bed, conversion from velocity to timescale in years (yrs) can be computed. Table 1 summarizes these results, which also are available in the associated data release (Burke, 2022). Case studies in the next three sections illustrate the application and interpretation of this quantitative migration modeling.

The influences of porosity, permeability, pressure, temperature, and saturation on migration velocity are investigated. According to this sensitivity analysis, the principle factors affecting migration velocities are the porosity and permeability of the carrier beds. Higher porosities and permeabilities correspond to increases in migration velocities and reductions in timescales. Pressure influences thermophysical properties of methane and CO₂; thus, a wide range of pressures from 0 MPa to 200 MPa (0–29,000 psi) are studied. Pressure also is accounted for in Darcy’s law. Higher pressures correspond to faster migration velocities and reduced timescales. However, changes in reservoir temperature along a carrier bed are expected to be scalar values and not exponential. Scalar changes in temperature would minimally affect density, viscosity, and P-wave velocity of the fluid, which would result in a negligible influence on migration velocities and timescales. Fluid compressibility asymptotically approaches a minimum value and remains invariant to temperature changes at these reservoir pressures. Investigation of various methane-water and CO₂-water fluid saturations are investigated at 25% increments, for example, 25% methane and 75% water. Results indicate that saturation influences fluid density, viscosity, and velocity in a linear fashion that is well below the resolution of this order-of-magnitude modeling technique.

The Haynesville Formation of the onshore U.S. Gulf Coast is an important domestic resource for petroleum. U.S. Geological Survey (USGS) assessments of undiscovered, technically recoverable resources estimate the Haynesville Formation contains means of 1.1 billion barrels of oil, a mean of 195.8 trillion cubic feet of natural gas, and a mean of 0.9 billion barrels of natural gas liquids (Paxton et al., 2018). Although the Haynesville is widely considered to be a self-sourced natural gas reservoir (Wang and Hammes, 2010; Hammes et al., 2011), the Haynesville and underlying Smackover source rocks are not thermally mature enough to generate the large volumes of dry gas produced along the southeastern margin of the Sabine uplift area in Louisiana (Paxton et al., 2021). However, Haynesville and Smackover Formation source rocks are generating dry gas in the adjacent North Louisiana Salt Basin (Torsch, 2012). Quantifying secondary migration velocities may reveal if the low-permeability Haynesville could represent a feasible migration pathway over a reasonable geologic timeframe.

Secondary migration velocities are modeled using average reservoir temperature (T = 177°C or 350°F) and representative rock properties (⁠$ϕ$ = 9%; $κ$ = 400 nD) of the Haynesville Formation (Table 2) after Wang et al. (2013). Accordingly, a reservoir temperature of 177°C (350°F) is used to condition pore fluid thermophysical properties over a broad range of reservoir pressures from 0 MPa to 200 MPa (0 to 29,000 psi). Modeling results (Figure 4d) indicate that interstitial pore fluid velocity is 3 × 10⁻¹⁰ m/s, which corresponds to 1 × 10⁵ years (100,000 years) for secondary migration to transit 1.0 km of these low-permeability carrier beds. The distance from thermally mature Jurassic source rocks in the middle of North Louisiana Salt Basin to the center of the productive trend of the Haynesville Formation in the Sabine uplift area is approximately 155 km (96 mi). Given this velocity and distance, 15.5 million years (Ma) are required for secondary migration of dry gas in this example pathway, which is beyond the production timeframe of transporting fluids from the reservoir into the wellbore but well within the geologic timeframe to charge a reservoir. Thus, the potential exists for additional hydrocarbon charges to have originated in an adjacent basin due to fluid flow based on Darcy’s law.

Geologic sequestration and long-term subsurface storage of supercritical CO₂ represent a solution for decreasing atmospheric emissions of this anthropogenic greenhouse gas (Benson and Cole, 2008). Identification of appropriate saline formations for CO₂ containment can be determined, in part, from the velocity and timeframe that fluids migrate through them. Supercritical CO₂ migration velocities are modeled (Burke, 2011) for a broad range of porosities and permeabilities likely encountered in the sedimentary sequence (Table 3). Based on the geothermal gradient and target formation depth, an average reservoir temperature of 93°C (200°F) is used (Burke, 2011) to constrain pore fluid thermophysical properties over reservoir pressures ranging from 0 MPa to 200 MPa (0–29,000 psi).

Modeling results (Figure 4) are shown as a function of porosity. For a storage formation with 1.0 D permeability (Figure 4a) and 10% porosity, results indicate that interstitial pore fluid velocity is 3 × 10⁻³ m/s, which corresponds to days to weeks for migration of these pore fluids to transit 1.0 km of the storage formation. Similarly, migration velocity through porous media with 1.0 mD permeability (Figure 4b) and 10% porosity is 3 × 10⁻⁶ m/s, which corresponds to 94.6 m/yr. This implies decades of containment for millidarcy permeability rocks exhibiting porosities approximately 10%. Migration velocity through 1.0 μD permeability (Figure 4c) strata at 10% porosity is 3 × 10⁻⁹ m/s, which corresponds to 9.46 × 10⁻² m/yr for containment spanning centuries to millennia. Migration velocity through 1.0 nD permeability strata (Figure 4d) with 10% porosity is 3 × 10⁻¹² m/s, which corresponds to 9.46 × 10⁻⁵ m/yr, representing several million years of containment; however, injections pressures required to sequester fluid in this strata are unfeasible with present engineering technology. Migration velocities in media with picodarcy permeabilities (3 × 10⁻¹⁵ m/s; not shown) could theoretically harbor fluids for billions of years, which exceed the tectonic cycle of most sedimentary basins. This migration model was developed to quantify supercritical CO₂ flow timescales (Burke, 2011) as a way to divide the subsurface into permeability-based flow units in support of the USGS methodology to assess the storage capacity of sequestered carbon dioxide into saline formations (Burruss et al., 2009; Brennan et al., 2010).

The Illinois Basin-Decatur Project is the world’s first large-scale carbon capture and storage demonstration project from a biofuel source in which 1.1 million tons of CO₂ are sequestered over a three-year period from 2011 to 2014 at an injection rate of 1000 tons per day (Gollakotaa and McDonald, 2014; Bauer et al., 2016). At this site, the Mesoproterozoic (1467 ± 25 Ma) basement is fractured rhyolite with a clay-rich matrix, and it is unconformably overlain by thin (approximately 33 m or approximately 100 ft) low permeability fine- to medium-grained sandstones (Leetaru and Freiburg, 2014), which are informally known as the Cambrian pre-Mount Simon interval (Freiburg et al., 2014). The 2130 m (7000 ft) deep target formation is the upper Cambrian Mount Simon sandstone, which is a 460 m (1500 ft) thick sand-rich fluvial system containing braid plain, flood plain, alluvium, and playa deposits (Leetaru and Freiburg, 2014). The regionally extensive middle to upper Cambrian Eau Claire Formation, comprised of low-permeability siltstones and mudstones, forms the 100–150 m (330–500 ft) thick primary seal preventing buoyant CO₂ migration into the overlying strata; the upper Ordovician Maquoketa Shale Group at 820 m (2700 ft) and the Devonian New Albany Shale at 650 m (2100 ft) form the shallower secondary and tertiary seals (Leetaru and Freiburg, 2014), respectively.

Four textural facies are interpreted within the Mount Simon, including unbleached fine-grained and coarse-grained sandstone beds, as well as bleached fine-grained and coarse-grained sandstone beds (Ritzi et al., 2016). Although the specific cause of bleaching in the Mount Simon is unknown (Ghose, 2017), Ritzi et al. (2016) observe permeability reduction and porosity loss due to pore cementation and mechanical compaction in the bleached intervals, as compared with unbleached intervals. Average porosity and permeability for each of these textural lithofacies are given in Table 4 (Ghose, 2017); the representative Mount Simon rock properties for our velocity modeling are $ϕ$ = 17% and $κ$ = 45 mD. Migration velocities are modeled through 1.0 km of isotropic porous media exhibiting these rock properties while using a reservoir temperature of 93°C (200°F) to condition the thermophysical properties of fluid density, fluid viscosity, fluid P-wave transit time, and fluid compressibility over reservoir pressures ranging from 0 MPa to 200 MPa (0–29,000 psi). The wellhead injection pressure is approximately 9.3 MPa (1350 psi) (Finley et al., 2013). A hydrostatic pressure of 9.8 kPa/m (0.433 psi/ft) in the reservoir at 2,130 m (7,000 ft) in depth yields a reservoir pressure of 20.9 MPa (3,031 psi) and a pressure differential away from the turbulent flow of the injection plume is 11.6 MPa (1,681 psi). This results in a migration velocity of 6.06 × 10⁻⁵ m/s, which equates to 1910 m/yr or a distance of 19 km (11.8 mi) in the 10 years after injection began; however, this upper bound estimate assumes no pressure decay. According to an analog carbon capture and storage project in Shenhua, China (Zhu et al., 2015), 20%–50% of this upper bound is more realistic and takes into account the geomechanical response of the formation after carbon sequestration concludes. Therefore, the CO₂ plume may have migrated 3.8 km (2.4 mi) and as far as 9.5 km (5.9 mi) away from the injection wellbore. USGS passive seismic monitoring of this site has shown no increase in the number or magnitude of microseismicity events either radially away from the wellbore or occurring at depths other than the injection zone (Kavern et al., 2015;Bauer et al., 2016). Monitoring for potential leakage includes inspection of wellhead equipment, pressure monitoring of the CO₂ plume, neutron tracer logging, and groundwater quality monitoring during the injection phase and extending 50 years after injection ends; the most likely leakage pathway is expected to be through aging wellhead and pipeline equipment and not a geologic pathway (Gollakotaa and McDonald, 2014).

The model examines wellbore, reservoir, and basin scales of investigation; nanoscale molecular interactions at the fluid-mineral interface are not taken into account. Navier-Stokes momentum equations are used to describe the fluid dynamics of the system; chemical dissolution into other fluid phases or chemical reactions at the fluid-solid interface are beyond the scope of this study. The constitutive laws rely on fluid properties as a function of pressure and a representative temperature. Large changes in pressure may influence fluid properties; thus, a wide range of pressures is accounted for in the model. Linear changes in fluid temperature, as may be expected from fluid migration originating deeper in the basin, would have a negligible influence on fluid density, fluid compressibility, and fluid viscosity, resulting in a negligible influence on the resultant migration velocities and timescales. The model assumes fluid transport through homogeneous media in which representative rock properties are used as inputs. This scenario does not specifically account for porosity heterogeneity or permeability variation within the carrier bed; variation in rock properties can be accounted for by taking the weighted average of the resultant velocities from a series of models each with relevant $ϕ$-$κ$ combinations. Note that fault and fracture networks can be approximated by higher porosities and permeabilities. These models examined the migration velocities of methane and CO₂, which have molecular masses of 16.04 and 44.01 atomic mass units (amu), respectively. Heavier hydrocarbons are expected to decrease the resultant migration velocity. For example, ethane and propane (30.07 amu and 44.1 amu, respectively) presumably plot between methane and CO₂ shown in Figure 4 velocity plots.

This study introduces a quantitative methodology for calculating migration velocities of supercritical pore fluids through porous media under reservoir temperature and pressure conditions. Constitutive laws and equations of state are used to obtain the thermophysical fluid properties used to condition a time-dependent derivation of Darcy’s law for calculating fluid transit through a wide range of porosity and permeability carrier beds likely encountered in the sedimentary sequence. Three case studies demonstrate the application of this migration velocity modeling. Fluid flow velocities of methane are modeled to determine that potential exists for additional hydrocarbon charge of the Haynesville Formation to have originated in and migrated from an adjacent basin. This migration model is further adapted for calculating velocities of supercritical CO₂ migration in porous media as a means to identify permeability-based containment formations appropriate for geologic sequestration and long-term storage of this greenhouse gas. The model also is used to estimate the 3.8 km (2.4 mi) distance that geologically sequestered CO₂ could be expected to migrate radially from the injection wellbore at the Illinois Basin-Decatur project, a large-scale carbon capture, and storage demonstration project. The model is applicable for a wide range of reservoir conditions (0–200 MPa or 0–29,000 psi) and rock properties (⁠$ϕ$ = 0%–95%; $κ$ = 10 D–0.1 nD); the model can be adapted to accommodate other temperature conditions and geologic scenarios.

The author acknowledges the journal for the opportunity to present this work and the U.S. Geological Survey (USGS) Central Energy Resources Science Center. Technical reviews of the manuscript were conducted by J. Birdwell, J. Herrick, J. Kavern, S. Paxton, J. Pitman, B. Shaffer, and P. Warwick. This research was funded by the USGS Energy Resources Program.

Data tables associated with this modeling include migration velocities and timescales for supercritical methane and carbon dioxide (Burke, 2022).

A biography and photograph of the author are not available.