A New Multiscale Model For Flow and Transport in Unconventional Shale Oil Reservoirs Tien Dung Le∗ and Marcio A. Murad Laborat´orio Nacional de Computac¸a˜ o Cient´ıfica LNCC/MCT Av Get´ulio Vargas 333, 25651–070 Petr´opolis, RJ, Brazil email:
[email protected] [email protected] July 26, 2018
Abstract A new multiscale computational model for two-phase gas/oil flow in a multiporosity shale formation composed of three levels of porosity associated with nano/micro pores and hydraulic fractures is rigorously constructed within the framework of the homogenization procedure in conjunction with the Discrete Fracture Modeling, where fractures are treated as (n-1) interfaces (n=2, 3) immersed in the domain occupied by the matrix. Effective equations are obtained by upscaling the microstructural information of the shale oil formation with matrix composed of three distinct solid phases: the organic matter, constituted by kerogen aggregates containing particles and nanopores with adsorbed gas, and the pyrobitumen network, also containing an organic solid with micropores filled by tight oil and dissolved gas, along with the inorganic solid composed of clay, quartz and calcite, assumed impermeable and free of adsorption. Such distinct solid phases are separated from each other by the network of interparticle pores. Together with the pyrobitumen such a network form the pathways for multiphase flow in the matrix whereas the kerogen aggregates are treated as disconnected inclusions playing the role of storage sites for adsorbed gas. The homogenization of the multiphase flow model of black oil type gives rise to new pressure and saturation equations with effective coefficients strongly correlated with the shale microstructure, volume fractions and total organic content (T OC). Constitutive laws for the effective hydraulic conductivity and retardation parameter, which captures adsorption of methane in the nanopores, are numerically reconstructed by solving the local cell problems arising from the homogenization procedure. In particular the partition coefficient is computed from adsorption isotherms rigorously constructed within the framework of the Thermodynamics of confined gases seated on the Density Functional Theory (DFT). The effective equations in the matrix resemble in form of a generalized black oil model coupled with the two-phase flow model posed in the subdomain occupied by the network of hydraulic fractures. A macroscopic model is obtained by averaging the mass conservation equation across the fracture aperture giving rise to reduced balance laws posed in a network of reduced (n1)-dimension (n=2, 3) supplemented by a source term arising from the jump in the oil/gas fluxes in the shale matrix. The resultant coupled Discrete Fracture/Multiscale model consists of a first attempt at constructing a rigorous correlation between the nature of the macroscopic multiphase flow equations and the local shalemicrostructure mainly characterized by the simultaneous presence of inorganic and organic matters with the latter containing nanopores. Numerical simulations of gas/oil withdrawal are performed to accurately predict hydrocarbon movement in stimulated shale oil formations.
1
Key words: Shale oil reservoirs, Pyrobitumen, Kerogen, Homogenization, Two-phase flow, Black oil formulation, Density Functional Theory, Gas Adsorption in nanopores, Hydraulic Fractures, Discrete Fracture Modeling
1
Introduction
For decades shales have been envisioned by the oil and gas industry as source rocks of hydrocarbons or barriers for their movement. More recently, with the rapidly increasing demand for global oil and gas resources, new challenges have emerged giving rise to the so-called exploitation of unconventional sourcereservoir systems such as tight-gas sands, coalbed methane, heavy oil, shale gas and shale oil [1]. Because of their large reserves such geological formations are becoming increasingly important and have appeared as alternative sources of hydrocarbon supplies worldwide. The terminology unconventional refers to the hydrocarbon-bearing formations that generally do not produce economic flow rates unless effectively stimulated. Techniques to enhance production such as advanced horizontal drilling and multi-stage hydraulic fracturing have enabled the oil industry to economically extract hydrocarbons from ultra-tight formations by creating hydraulic fracture networks that connect huge reservoir surface area to the wellbore [2]. Shale-oil formations are characterized by microstructure and geochemical properties containing high degree of compositional heterogeneity, extremely tiny pores of nanometer size and ultra low permeability [3, 4, 5]. The complexity underlying the multiscale description of such geological formations lies in the presence of a nano-scale pore network within the organic matter, the interparticle pores which separate the inorganic and organic matters and the presence of natural and hydraulic fractures. In particular, the organic solid is composed of two different substructures commonly referred to as kerogen aggregates containing gas-wet nanopore network filled by adsorbed gas and pyrobitumen composed of post-oil solid bitumen and micropore network capable of hosting tight oil with dissolved gas and free gas [6, 7, 8, 9]. Experimental evidences have suggested that effective multiphase flow in the shale matrix occurs in the networks of interparticle pores and micro sized-pores in the post-oil solid bitumen whereas the gas-wet nanopores within the kerogen aggregates play the role of storage sites for methane adsorption due to the non-connected kerogen network often quantified by low values of the total organic carbon (T OC < 10%) [6]. In Fig. 1.1, we depict the natural distinct length scales addressed in the multiscale model proposed herein. As quoted before the stimulated region of the shale-oil exhibits four distinct types of pores along with four separate length scales. At the nanoscale, gas adsorption occurs in the nanopores of the kerogen. At the microscale the system is treated as a three-phase medium composed of impermeable inorganic matter, organic clusters (with secondary levels of porosity associated with the kerogen aggregates and the post-oil solid bitumen) and finally the network of interparticle pores. The aforementioned issues regarding shale morphology are of utmost important and need to be precisely incorporated in macroscopic hydrodynamic models. Accurate relations between the theories at different length scales have been previously addressed by the authors [10, 11] for single phase flow in shale gas reservoirs where effective nonlinear parabolic equations were constructed to describe coupled flow in the matrix and fracture subdomains. In this context, the inclusion of anomalous behavior of the sorbed gas in the nanopores of kerogen may be treated in a simplified manner invoking the classical Langmuir isotherm, whose accuracy is restricted to the monolayer adsorption picture when the characteristic length of the nanopores l is much larger than the diameter of the methane molecule d [10]. In the case of ultra nanopores, where for instance l < 5d, more sophisticated theories such as molecular dynamics, Monte Carlo simulations and Density Functional Theory (DFT) based on the thermodynamics with potentials exhibiting functional dependence on the local density profiles. The DFT framework has been previously adopted by the authors to compute the partition coefficient associated with methane adsorption in the nanopores [10, 11].
2
The development of accurate computational models for multiphase flow in ultra low permeability shale oil reservoirs still poses significative challenge issues. Most of the existing models are seated on classical black oil and compositional-type models with the magnitude of the input coefficients adapted to the low permeability and porosity scenarios and to accommodate sorbed gas via Langmuir adsorption isotherms and partition coefficients [12]. Moreover no representation in the flow equations has been accomplished so far to explicitly capture the distinction between flow and adsorption in the kerogen and post-oil solid bitumen. Despite the substantial improvement achieved in the development of multiscale models for single phase flow in shale gas, the extension of the methodology to multiphase flow remains an open question. In this work we provide the first steps toward the development of a sophisticated multiscale model for multiphase flows in shale oil reservoirs. To accomplish this task we build-up the macroscopic governing equations by homogenizing the microstructural local phenomena which take place in the complex local networks, whose composition and pore-size distributions are somewhat related to the history of the diagenetic processes [6, 7, 3, 4, 9]. We develop a new generalized black oil type-model for two-phase oil and gas flow in shale oil reservoirs by rigorously upscaling the local behavior of the two-phase flow model. Our microscopic picture considers the networks of interparticle pores separating inorganic and organic matters with the latter constituted by fracture-filling post-oil solid bitumen containing micro-pores saturated by tight oil along with disconnected kerogen inclusions playing the role of storage sites for gas adsorption. Based on the aforementioned methodology our main goal relies upon the construction of precise multiscale representations capable of correlating the macroscopic response of the medium with the complex behavior of the microstructure. In a similar to [10, 11], we proceed within the context of the DFT to describe the partition coefficient governing adsorption of the anomalous gas. The constitutive dependence of this parameter is then coupled through interface conditions with the movement of the oil/gas mixture in the interparticle pores and postoil bitumen, governed by the classical black oil model. The coupled black oil-type model for flow in the networks of interparticle pores and post-oil solid bitumen is upscaled to the homogenized (meso) scale, where the shale matrix is treated as a single continuum with mean flow and transport ruled by effective equations. The formal technique adopted for describing phenomena at the different length scales is based on the formal homogenization procedure seated on matched asymptotic expansions [13]. The numerical solution of the local microscopic cell problems allows to construct numerically the new constitutive laws of the effective coefficients (retardation coefficient and hydraulic conductivity). In particular the computation of the effective conductivity requires knowledge of the tortuosity function whose estimates can be obtained by adopting the self consistent homogenization method [13]. After establishing the multiphase flow model in the matrix we postulate governing equations in the hydraulic fractures and proceed by averaging the conservation laws across the fracture aperture. Thus, homogenized governing equations in the matrix are coupled with the mean flow posed in a reduced domain of (n-1)-dimension (n=2, 3) occupied by the fractures giving rise to a new expanded coupled problem which can be embedded in the framework of Discrete Fracture Modeling (DFM) posed simultaneously in the fracture/matrix domains [14, 15]. The subsystems in each subdomain are discretized by the finite difference method and the nonlinear coupling between them is treated iteratively by fixed point Picard-type schemes. In contrast to the pressure equation, the saturation problem is solved explicitly in time, with time step constrained by the CFL condition, and numerical fluxes computed by the first-order upwind scheme. Numerical simulations of gas/oil withdrawal in a stimulated shale oil reservoir are presented for the particular arrangement of parallel fractures orthogonal to a production horizontal well. We adopt a structured mesh in the fracture discretization with reduced interface elements coupled with rectangular cells in the partition of the shale matrix domain. In this scenario our set of numerical simulations illustrates the influence of parameters such as fracture/matrix permeabilities, pyrobitumen content, T OC, and porosities upon pressure, flux, saturation profiles and production curves. The accuracy of our computational simulations relies on the multiscale microstructure-based nature of the model, which is capable of providing a detailed sensitivity analysis of the response of unknown 3
with respect to perturbation in the input parameters. -3
O(1m)
O(10m)
PHydraulic PPFracture
Matrix
Ω*F Matrix Matrix
Matrix
ΩM
ΩI PNanoporesP PofPkerogen
Pyrobitumen
Kerogen
-9
Ωp
ΓmI
Ωm Γmp
Γ
Ωkmk
ΩI
InterparticleP PPPPPpore
InorganicP PPmatter
-6
O(10m)
O(10m)
Figure 1.1: Microstructural configuration associated with nano, micro, meso and macroscales in a shale oil formation
2
Microscopic model
At the microscale, the shale-oil is composed of three distinct solid phases; the kerogen aggregates, which are envisioned as a porous medium with local organic particles and nanopores filled by adsorbed gas, along with the pyrobitumen filled by tight oil with dissolved gas [8] and the inorganic solid assumed impermeable and free of adsorption. Three solid phases are separated from each other by the network of micrometer sized pores referred herein to as the system of interparticle pores [10]. On the other hand, following [8], besides of acting as storage sites for tight oil, the pyrobitumen network along with the interparticle pores provide pathways for the mean flow in the matrix. In our subsequent developments, we already start with one level of averaging in order to overcome the tremendous complexities associated with upscaling the Stokes problem in two-phase flow with unstable interfaces. Thus we begin by postulating mean equations strictly taken across the pore space filled by the fluid mixture with the exception of the adsorbed gas in the kerogen. The fluids in the interparticle pores and pyrobitumen are in the bulk state and therefore we begin by postulating the blackoil formulation for the two-phase flow in the domains occupied by the interparticle and micrometer sized pores of the pyrobitumen. For simplicity and without loss of generality gravitational effects are neglected. Let ΩI , Ωp and Ωk ⊂
