CO₂, as one of the greenhouse gases, plays a key role in the cause of global warming, and its accumulation in the atmosphere is still increasing (Allen et al. 2014). As a result, CO₂ mitigation has attracted increasing attention (Edenhofer et al. 2014). CO₂ capture and sequestration (CCS) technology is an alternative to effectively reduce the CO₂ emission in the short term (IEA 2013). Popularly, brine or saline aquifers, oceans, depleted oil reservoirs, and coal beds are considered potential geological formations for CO₂ storage. For example, Fujii et al. (2010) and Zahid et al. (2011) argue that saline aquifers are the most viable alternative to store CO₂ because they possess the greatest potential for carbon storage and wide geographical spread. Other studies suggest that CO₂ is a working fluid that can be used to enhance oil recovery (EOR) and natural gas recovery (EGR) (Hughes et al. 2012; Liu et al. 2015a, b; Middleton et al. 2015). Injecting CO₂ into oil or gas reservoirs is promising because it can also offset the costs of CCS (Koide et al. 1992; Blunt et al. 1993; Gunter et al. 1997). When CO₂ is injected into a deep geological formation in the liquid or supercritical state, it will cause large volumes of formation fluids such as oil, water (brine) or natural gas to be physically displaced (Kazemifar et al. 2016; Wang et al. 2013; Zhang et al. 2011b). In the CO₂ displacement process, the major concern is the primary CO₂ plume migration, which is closely related to how much CO₂ can be stored in the respective porous subsurface sedimentary formation, the displacement efficiency of the oil/gas and the security of the stored CO₂ (Berg and Ott 2012; Chen and Zhang 2010; Song et al. 2014). Therefore, a better understanding of the displacement mechanisms of CO₂ and formation fluids in porous media is essential to assess the CO₂ leakage risks and predict the amount of CO₂ that can be absorbed and the oil/gas production as well.

Predicting changes in the flow and transport properties of CO₂ and formation fluids in porous rock is still a challenging issue owing to the multi-phase flow, complexity of the fluid–rock interactions and intrinsic heterogeneity of porous rocks (Garcia-Rios et al. 2015). Thus, the CO₂ displacement process in porous media has attracted increasing attention. Extensive research on CO₂ displacement has been performed over multiple spatial scales that range from a few nanometres to tens or hundreds of kilometres, and it is also clear that the appropriate methods to be employed vary between different scales. At the macro-scale, average properties must be employed, as an explicit description of the pore-scale phenomena is simply impossible. Additionally, the physical and chemical processes occur in discontinuous pore geometries and are controlled by interfacial processes, which have a large influence on the large-scale phenomena of the system and ultimately require a pore-scale perspective for obtaining a better understanding of the mechanisms (Blunt et al. 2013; Wildenschild and Sheppard 2013; Ferrari 2014; Morais et al. 2016). In the last decades, the pore-scale modelling of CO₂-formation fluid displacement has been promoted by three factors: the recent advances in porous media characterizations by visualization techniques, some novel numerical theories, and advances in high-performance computing (Ferrari 2014). However, studies on the CO₂-formation fluid displacement using these promising technical and theoretical approaches have not been systematically analysed in the published literature.

Thus, the objective of this contribution is to outline the recent developments and advances in the simulation and visualization of the CO₂-formation fluid displacement process occurring at the micro-scale and nano-scale—also referred to here as the pore-scale—and to highlight perspectives for future research. In Sect. 2, some relevant important processes and characterizations of the CO₂-formation fluid displacement are briefly discussed. In Sect. 3, several novel visualization techniques and numerical simulation methods for the CO₂-formation fluid displacement are introduced in detail, including a discussion of their strengths and weaknesses. Finally, Sect. 4 provides a comprehensive overview of the applications of CO₂ displacement to enhanced shale gas production including the status of research, challenges and possible solutions.

2.1 Viscous fingering

In porous media, viscous and capillary forces play key roles in impacting the CO₂-formation fluid displacement. The injected CO₂ viscosity is generally lower than that of formation fluids (Nordbotten et al. 2005), and the interfacial tension between the CO₂ and formation fluids depends on the pressure, temperature, and other system conditions (Chiquet et al. 2007; Espinoza and Santamarina 2010). Two dimensionless parameters, namely, the capillary number (Ca) and viscosity ratio (M), are commonly employed to describe the various forces encountered by the two fluids during their displacement in the porous media. These two parameters are defined as

where \(\mu_{1}\) and \(\mu_{2}\) are the viscosity of the displacing fluid and the viscosity of the displaced fluid, respectively. \(V_{2}\) is the bulk velocity of the displacing fluid, \(\gamma\) is the interfacial tension between the two fluids, and \(\theta\) is the fluid–fluid contact angle with the solid surface.

Based on the values of Ca and M, three displacement patterns can be identified, as shown in Fig. 1: viscous fingering, capillary fingering and stable displacement (Lenormand et al. 1988; Zhang et al. 2011a; Yamabe et al. 2015). In the CO₂ displacement process, the fluid flow speed is low in subsurface rocks (Ca ≪ 1), and the viscosity ratio is considered to be less than 1 (M < 1) (Yamabe et al. 2015). Hence viscous fingering can commonly occur during the CO₂-formation fluid displacement. Viscous fingering is the formation of special fingers with morphologically unstable interfaces of the invading fluid. A deep understanding of viscous fingering is essential for accurately predicting the migration and transport of the injected CO₂ and formation fluids within the pore structure of the formations.

Classification of displacement patterns based on Ca and M. The boundaries marked with blue lines are those noted by Lenormand et al. (1988), whereas those marked with red lines are those from Zhang et al. (2011a). Modified after Yamabe et al. (2015)

2.2 CO₂–rock interaction

With a deep understanding of the CO₂-formation fluid displacement process, the CO₂–rock interactions intervene in more areas than initially anticipated, whereby certain areas have attracted less attention to date. When CO₂ is injected into formations, scCO₂ will dissolve, and interaction will occur between the CO₂ and formation rock (Gaus 2010). As a consequence, these interactions will lead to changes in the pore structure of the formation rock, decreasing the porosity and permeability of the rock. These changes in the pore geometry of the formation rock will seriously impact the fluid transport in the CO₂-formation fluid displacement. Thus, there are three key issues summarized in the review of Gaus (2010), that is, where do CO₂–rock interactions occur, what are the drivers of CO₂–rock interactions and what is their influence on the porosity and permeability, and how can the CO₂–rock interactions be assessed.

3.1 Experimental visualization

Laboratory experiments on CO₂-formation fluid displacement can provide valuable insights into processes that control the displacement process at the pore scale of porous media and improve our understanding of multi-phase flow. Pore-scale experiments also permit the evaluation of the constitutive relationships that are used in large-scale simulators (Dehoff et al. 2012). Roughly, there are three main visualization approaches to detect and observe CO₂ displacement processes at the pore scale: magnetic resonance imaging (MRI), X-ray computed tomography (CT) (Shi et al. 2009, 2011a, b; Zhao et al. 2011a; Alemu et al. 2013; Berg et al. 2013b) and fabricated micromodels (Ferer et al. 2004; Cottin et al. 2010; Zhang et al. 2011a; Wang et al. 2013; Al-Housseiny et al. 2012). These methods are recognized as effective approaches to directly characterize the transport properties of the displaced fluid and displacing fluid at the pore scale.

3.1.1 X-ray CT

X-ray CT is a non-invasive and non-destructive imaging technique used to characterize the internal physical structure of the sample and allow for the three-dimensional (3D) in situ visualization of the fluid phases (Lindquist et al. 1996; Hunt and Bajsarowicz 1988; Wellington and Vinegar 1985, 1987; Schlüter et al. 2014). The theories and devices of X-ray CT are described by Wildenschild and Sheppard (2013) in details. In the last decades, X-ray CT has been widely used to investigate CO₂-formation fluid displacement. For instance, Berg et al. (2010) investigated the miscible displacement performance of a CS₂–nC₁₀ system at 1.5 MPa and 25 °C using CT scanning. The authors concluded that the front fingering could affect the mixing zone area. Berg et al. (2013a) also studied the displacement and mass transfer between saturated and unsaturated CO₂-brine systems in sandstone using CT, and 3D saturation patterns during the injection were obtained (Fig. 2). Seo (2004) and Seo and Mamora (2005) conducted the first experimental study of enhanced gas recovery by injecting CO₂ through a dry carbonate core that was saturated with CH₄ at 20–80 °C and 3.55–20.79 MPa. The distribution of CH₄ and CO₂ in the core was obtained using an X-ray CT scanner at the breakthrough time. The recovery of CH₄ was 73%–87%, and the dispersion coefficient of CO₂ was approximately 0.0–0.12 cm²/min. In addition to the above studies using real reservoir rock cores, some researchers chose sand packs as an experimental material for basic research on the displacement process. Liu et al. (2015b) investigated the CO₂–CH₄ displacement process in sand packs, and the displacement process was scanned with a resolution of approximately 34.2 mm by an X-ray CT scanner. Concurrently, the dispersion coefficients were obtained, and the effects of the temperature, pressure, CO₂ injection flow rate and mean diameter of the glass beads were analysed. Through the X-ray CT scan, the detailed displacement process in the sand pack sample was obtained noninvasively, as shown in Fig. 3. The existence of the mixing zone was proven, and no interface appeared in the process of CO₂–CH₄ displacement.

3D saturation patterns (CO₂ in red/yellow), a unsaturated CO₂-brine systems and b saturated CO₂-brine systems. The rock matrix is visible as a semi-transparent background. The displacement shows the gravity over-run of CO₂. (Berg et al. 2013a)

a CT scan images of the same vertical section in displacement in the BZ04 glass bead sand pack sample with time; b corresponding CT scan images of the same middle cross-section located on the red line of the vertical section of (a); The colour scale is set based on grey scale, with glass beads painted as white with edges that may be red or yellow, blue and dark blue are CH₄, green represents CO₂, and cyan indicates a mixture of CH₄ and CO₂. (For the interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) (Liu et al. 2015b)

However, based on previous studies, the front fingering and the mixing zone of the displacement process can only be qualitatively analysed by the CT imaging and cannot be quantified. Moreover, the microscopic fluid state cannot be detected using medical CT scanners owing to the insufficiently high resolution (millimetre-scale) compared with the pore sizes, which could make the results of the experiments less credible. The resolution of nano-CT (50 nm)is high, but it takes a long time to collect one image unless fast synchrotron-based sources are used.

3.1.2 MRI

MRI is similar to X-ray CT in that non-invasive visualizations of the fluid movement and distribution can be obtained in the core. The differences are that unprecedented quantitative information about the distribution of fluid sand information about the rock structure corresponding to local regions within the porous media can be provided by MRI during the displacement processes (Song et al. 2012). Special features of MRI include its abilities to measure the fluid flow velocities and to distinguish between different liquids by the differences in their intrinsic MRI properties, which are difficult to obtain by other methods. It also has a shorter measurement time and is much less hazardous in comparison with X-ray CT (Xiong et al. 2016). Using the MRI technique, some studies were performed on the CO₂-formation fluid displacement in porous media. Suekane et al. (2005, 2006, 2009) used the MRI technique to directly visualize the distribution of supercritical CO₂ (scCO₂) in a bead-pack core containing water at a certain temperature and pressure that approximate aquifers at a depth of approximately 1000 m. Liu et al. (2011) analysed the immiscible and miscible CO2 displacement of oil in glass bead packs using MRI technology. In their study, the local Darcy velocities of CO₂ were acquired based on the oil saturation data obtained from MRI images, and the effects of the capillary pressure, viscosity and buoyancy on the CO₂ displacement process can also be evaluated using a special core flood analysis. Zhao et al. (2011a, b) proposed a method to calculate the average velocity of the CO₂ front and implemented CO₂–oil miscible displacement in a sand pack using MRI. The local Darcy phase velocities of the CO₂ and oil were calculated with this method and are shown in Fig. 4. Song et al. (2012) studied the CO₂-water displacement process in packed glass beads using high-resolution MRI. MRI images of the CO₂ displacement were obtained, and the effects of the displacement type, CO₂ injection flow rate and effective porosity were considered. A summary of the results of the displacement experiments is shown in Table 1. Song et al. (2014) also analysed the in situ mixing zone performance of CO₂-oil miscible displacement flows in a sand pack using MRI. The mixing zone length, CO₂ frontal velocities, longitudinal dispersion coefficient and Peclet number were quantified based on magnetic resonance imaging. Yang et al. (2015) investigated the dynamic stability characteristics of fluid flow in CO₂-oil miscible displacement using an MRI apparatus, and the CO₂ frontal velocities and mixing zone length were visually quantified.

Local phase velocities of a CO₂ and b oil. (Zhao et al. 2011b)

Although a number of studies have been conducted to investigate the CO₂-formation fluid displacement in a sand pack using MRI, and some CO₂-formation fluid displacement mechanisms have been clearly analyzed, there have been almost no studies on natural rocks owing to the harsh conditions of MRI against the experimental materials. As the geometry of the natural rock sample is more complex than that of the sand packs, the reliability achieved for the previous experimental results can be compromised.

3.1.3 Fabricated micromodels

The rock cores have an advantage in characterizing individual formations, but on the pore scale, the fluid flow is difficult to monitor since sophisticated and unique micro-tomographic facilities are needed to visualize the internal distribution of the fluids within the rock cores. The natural media also presents other challenges in that the porosity, pore size, connectivity, and wetting properties are unable to be independently manipulated. These limitations can be overcome by micromodels, which allow for the visualization of the fluid distribution using cameras with or without fluorescence microscopy. The subsequent quantification of the fluid transport and interfacial area may provide mechanistic insight into the physical fluid displacement process at the microscopic level. Experimental studies on the CO₂-formation fluid displacement have been widely conducted in 2D micromodels to reveal the mechanisms of the fluid displacement.

Er et al. (2010) investigated the porous matrix and fracture interaction during CO₂ injection into a glass micromodel initially saturated with oil, demonstrating the importance of the CO₂–oil interaction near the matrix-fracture interface. Chalbaud et al. (2009) visualized the scCO₂ displacement of water in glass micromodels, but the saturations of the fluids could not be quantified owing to difficulties in distinguishing between the fluids. Riazi et al. (2011) simulated CO₂-oil and CO₂-water displacement processes in an etched glass micromodel, and the results showed a faster CO₂ breakthrough in the micromodel saturated with oil than that with water. The above studies were all limited to qualitative visualizations of the fluid flow during the displacement process, and quantitative information about the interfacial areas and fluid saturations was not acquired. The study of Zhang et al. (2011b) was the first to quantitatively evaluate the influence of the capillary forces and porous media heterogeneity on the liquid CO₂-water displacement in a dual-permeability pore network micromodel using fluorescence microscopy. Wang et al. (2013) simulated the effects of the injection rates and injection methods on the scCO2 displacement of water in a homogeneous micromodel, and the results of the simulation are quantitatively compared using the obtained fluorescent images. Kazemifar et al. (2016) quantified the flow dynamics of the scCO2 displacement of water in a 2D porous micromodel by combining fluorescence microscopy and microscopic particle image velocimetry (PIV).

Unfortunately, micromodels lack the complex geometry of real media, which often have multiscale and random characteristics that will dictate the fluid and solute transport. In addition, these experimental studies used simple two-dimensional porous media. To overcome the above limitations, Ju et al. (2014a) incorporated several advanced technologies, such as CT scan, three-dimensional (3D) reconstruction (Gomi et al. 2007; Geiger et al. 2009; Hajizadeh et al. 2011; Zhang et al. 2013; Birk et al. 2014; Ju et al. 2014b), and 3D printing, to produce a physical model representing the natural rock. A 3D digital image can be obtained to represent the real rock structure based on the high-resolution X-ray CT imaging and the reconstruction, and then its three-dimensional micromodel using transparent material can be generated using a 3D printer. This novel idea offers a very good opportunity to accurately investigate the properties of complex flows.

These experimental simulation methods greatly promote our understanding of the CO₂-formation fluid displacement mechanism. However, it is worth noting that these simulation experiments are performed based on the pore geometric characteristics of a convention formation that has a minimum pore size on the micrometre scale. The pore sizes and pore throat radii of unconventional gas reservoirs range from 1 to 300 nm, which are much smaller than those of conventional reservoirs with pore sizes in the range of 1–100 μm (Cipolla et al. 2009). These popular experimental simulation methods are thus not able to investigate fluid flow at the nano-metre scale because the maximum precision of these techniques is on the micro-meter scale in general. Special characteristic features that can be produced at the nanometre scale include a high capillary pressure, low porosity, and high wetting phase residue saturation (Wu 2014). Other transportation mechanisms may also occur. The displacements of two-phase flow in nano-pores are poorly understood. Wu (2014) developed a lab-on-chip approach for the direct visualization of the fluid flow behaviour in nanoscale channels using an advanced epi-fluorescence microscopy method combined with a nanofluidic chip. This method is expected to improve the understanding of the CO₂ displacement behaviour on the nano-scale. However, it is still a qualitative analysis method and cannot quantifiably reflect the displacement process.

3.2 Numerical simulation

A large number of the porous media studies on fluid displacement have been numerically performed, and numerical simulations can provide an economical and efficient pathway to explore the influences of the flow and physical parameters in various complicated porous media. Pore-scale simulations provide a level of information on flow characteristics that cannot be obtained in laboratory experiments, for example, the full pressure and velocity fields and the position of the interface. Modelling approaches of multiphase flow at the pore scale can be divided into two categories: direct simulations, in which the flow equations are directly solved on a discretized pore space obtained from images of rock cores, and network modelling, in which simplified flow equations are solved in an idealized pore network extracted from the real geometry (Blunt et al. 2013).

Since the first pioneering work (Fatt 1956), pore network models have been widely used to simulate multi-phase fluid flow at the pore scale. However, the validity of the pore network models is always limited by the approximations and simplifications used (Ferrari and Lunati 2013). Many pore network models may not be suited for the simulation of multiphase fluid flow in complex porous media (Crawshaw and Boek 2013). One of the most popular approaches is the lattice Boltzmann method (LBM), in which particles move and collide on a discrete lattice in such a way that the average motion of the particles mimics the solution of the Navier–Stokes equations (Shan and Chen 1993; He and Luo 1997). Since the late 1980s, the LBM has emerged as a powerful tool for numerical simulations and investigations of a broad class of complex flows. The LBM can be considered as a mesoscopic method, occupying the middle ground between the microscopic molecular dynamics and macroscopic fluid dynamics. Unlike pore-network models, which use a simplified representation of the pore geometry and approximate transient flow with a steady-state Poiseuille law, LBM can simulate complex multiphase flows in pore media with natural pore geometries (Bandara et al. 2013). The basic idea of the LBM is to construct simplified kinetic models that incorporate the essential physics of microscopic or mesoscopic processes to ensure that the macroscopic averaged properties obey the desired macroscopic equations (Ghassemi and Pak 2011). LBM has many advantages over the conventional grid-based computational fluid dynamics (CFD) methods, such as volume-of-fluid (VOF) and level-set (LS) methods, especially in handling complex boundaries, flexible reproducing interfaces, incorporating the microscopic interactions between multiple phases, and the parallelisation of the algorithm (Liu et al. 2015a).

Based on the above analysis, at present, the numerical simulation of the CO₂ displacement in porous media is mainly performed using the LBM, and few studies use conventional methods to simulate the displacement process. Gunde et al. (2010) investigated the displacement of oil by CO₂ within the pore structures of a realistic porous geometry (730 μm × 450 μm) derived from processed micro-CT images based on the finite element method. However, the COMSOL simulations require a minimum of 48 h of CPU time to obtain a converged solution for a 10 s displacement process. Berg and Ott (2012) used the phenomenological two-phase extension of Darcy’s law to study the stability of CO₂-brine immiscible displacement in a two/three-dimensional numerical grid. They found that the capillarity, viscosity ratio, gravity and length scales have effects on the stability of the displacement.

Liu et al. (2014) employed LBM to simulate liquid CO₂ (LCO₂) displacing water in a dual-permeability pore network (Fig. 5) and randomly heterogeneous pore network (Fig. 6). The extent and behaviour of the LCO₂ preferential flow (i.e., fingering) is found to depend on the capillary number (Ca). In comparison to the dual-permeability case, in the heterogeneous pore network case, the transition from capillary fingering to viscous fingering occurs at a higher Ca, and the LCO₂ saturation is higher at a low Ca but lower at a high Ca. Yamabe et al. (2015) simulated the supercritical CO₂–water drainage displacement using LBM in a 3D rock model (1 mm × 1 mm × 1 mm) that was constructed by a random packing process with a number of spherical grains at a fixed viscosity ratio and at various capillary numbers (by changing the injection pressure, which induces changes in the flow velocity). It was observed that the viscous fingering was dominant at high injection pressures, whereas crossovers of viscous and capillary fingerings were observed, accompanied by Haines jumps, at low injection pressures. The capillary number and displacement mechanisms are dependent on the geometry of the porous media. Then, they (Yamabe et al. 2016) further evaluated the relationship of V _P–CO₂ saturation, which is the most important relationship in the geophysical monitoring of CCS projects, applying dynamic wave simulation with the FD approach to LB results and the 3D rock model. A lower V _P was observed in the high capillary number case than at a low capillary number at the same saturation. The difference in the V _P response to CO₂ saturation is caused by the CO₂ distribution features. Yang et al. (2015) also used a three-dimensional LBM to simulate the pore-scale scCO₂-oil miscible displacement in upward vertical displacements. At a low kinematic viscosity ratio, a front fingering was found, and a CO₂ breakthrough occurred quickly (Fig. 7a). At a large kinematic viscosity ratio, the displacement was stable with a piston front. The high kinematic viscous oil restrained the buoyancy of the supercritical CO2 and impeded the displacement with a pore-scale back flow, which might lead to a low oil recovery factor (Fig. 7b). Their study also indicated that the agravic miscible displacement is slightly quicker than the miscible displacement with gravity in upward miscible displacements.

Final fluid distributions in the dual-permeability microfluidic system at various capillary numbers: a log Ca = −4.36, b log Ca = −4.16, c log Ca = −4.06, d log Ca = −3.59, e log Ca = −3.36, f log Ca = −3.06. LCO₂ and water are shown in red and blue, respectively. (Liu et al. 2014)

Final fluid distributions in the randomly heterogeneous system at a log Ca = −4.36, b log Ca = −4.16, c log Ca = −4.06, d log Ca = −3.59, e log Ca = −3.36, f log Ca = −3.06. LCO₂ and water are shown in red and blue, respectively. (Liu et al. 2014)

Simulation results in miscible displacement flows. Kinematic viscosity ratio of CO₂ to oil is a 1:1 and b 1:15 (Yang et al. 2015)

From the above review, we see that only a limited number of numerical simulation studies have been reported focusing on fundamental processes at the microscopic pore scale during CO₂ injection into a porous formation. These studies also did not systematically investigate the CO₂ displacement mechanisms. During the CO₂ flow in saline aquifers or oil/gas reservoirs, the main physical and chemical interactions among the water, oil/gas and rocks include (Ju et al. 2013): (1) the immiscible and miscible displacement of oil, water and CO₂; (2) convection and diffusion; (3) the phase change and behaviour of insitu fluids in flow and mass transfer; (4) precipitation/dissolution and its effects on porosity and permeability; (5) the CO₂ adsorption characteristics in the adsorption of the reservoir. Therefore, in further work, these interactions need to be considered. In addition, although there have been a few studies on the use of LBM for the simulation of the CO₂ displacement process in 3D porous media, the 3D pore structure of real rock has not been used. The main reason is that attaining an accurate description of the 3D microstructure of porous rock is difficult, laborious, time-consuming, and economically expensive with the use of advanced 3D CT measurement techniques, e.g., synchrotron X-ray CT (Coker et al. 1996; Spanne et al. 1994), and directly using the 3D microstructure of porous rock for the numerical simulations will entail a massive computational cost. Recently, Ju et al. (2014b) proposed a faster, more accurate, and more efficient reconstruction algorithm (an improved simulated annealing algorithm reconstruction method) to create a 3D model of the porous rock microstructure. On the other hand, the appearance and application of modern high-performance computers offer a new way to study complex flows in a domain with realistic pore geometries.

Shale gas, as a renewable energy source, has become more and more important to the worldwide energy supply owing to its low pollution, abundant reserves, and wide distribution (Curtis 2002; Burnham et al. 2012). It usually has two states stored in shale, that is, free gas in natural cracks and adsorbed gas on the organic matter surfaces (Yin et al. 2016). Generally, a significant part of the shale pores are distributed on the nanometre scale, which results in a large internal surface area (Loucks et al. 2009; Kuila and Prasad, 2013; Mosher et al. 2013). Thus, a large proportion of the shale gas is in an adsorbed state in the nanopores (Chareonsuppanimit et al. 2012). A large number of nanometre pores also lead to the extremely low permeability of the shale (Javadpour et al. 2007; Javadpour 2009; Heller et al. 2014; Hinai et al. 2014). There are many difficulties in shale gas extraction based on the above two issues, and effective methods are urgently required to enhance shale gas recovery (ESGR). Hydraulic fracturing is currently the only fracturing technology utilized in commercial shale gas extraction (Yin et al. 2016). However, water-based fracturing fluids have many drawbacks, including aggravating water scarcity and environmental problems (King 2010; Nagel et al. 2013; Rutqvist et al. 2015). Recently, scCO₂ has attracted widespread attention as a fracturing fluid for use in ESGR (Rogala et al. 2014). ScCO₂ has a lower viscosity and similar density to water (Ishida et al. 2012) and its surface tension is almost zero, so in the nanometre pore structure of shale, it could easily flow and induce more complicated fractures. CO₂ injection has also proven to be a potential way to displace shale gas in shale gas formations because the affinity of CO₂ for the shale surface is stronger than that of shale gas (Kang et al. 2011; Heller and Zoback 2014). In addition, shale gas reservoirs are both reservoirs and cap rock, which can realize the effective storage of carbon dioxide. Therefore, the integration of scCO₂-enhanced shale gas production and geological sequestration is proposed (Fig. 8). When scCO2 is injected into shale gas reservoirs, four integration effects are generated and three scientific questions need to be solved. A number of studies on the CO₂ fracturing mechanism have been reported (Ishida et al. 2012; Li et al. 2015; Huang et al. 2015; Bennour et al. 2015; Chen et al. 2015), and this paper does not provide a detailed introduction. However, the latter two issues are related to the displacement of multiphase fluid, and there are few studies in these areas owing to the special properties of shale gas reservoirs. The existing problems and possible solutions for studying the CO₂ displacement mechanism in the shale pore structure will be introduced in detail below.

Integration of scCO₂-enhanced shale gas production and geological sequestration

4.1 Special properties of shale reservoirs

4.1.1 Nanopore structure

Unlike conventional reservoirs, shale reservoirs are usually dominated by nanometre-scale pores (Nelson 2009; Loucks et al. 2009; Kuila and Prasad 2013; Mosher et al. 2013). According to a previous review on the simulation of CO₂-formation fluid displacement, at present numerical simulation is the only feasible method to investigate the CO₂ displacement process in nanometre-scale pores, and the quantitative characterization of the shale nanopore structure and pore network modelling are required before the numerical simulation of the CO₂-shale gas displacement can be performed. For this end, several imaging techniques and modelling methods are utilized to characterize the shale pores, such as nano-CT, field emission and scanning electron microscopy (FE-SEM), transmission electron microscopy (TEM), and focused ion beams and scanning electron microscopy (FIB-SEM) (Loucks et al. 2009; Bernard et al. 2012; Chalmers et al. 2012a, b; Curtis et al. 2012; Milliken et al. 2013). Of these methods, the maximum spatial resolution of nano-CT is only 50 nm, and so it cannot satisfy the requirements for shale nanopores of less than 50 nm, which can make the prediction of the fluid migration properties less credible. The spatial resolution of TEM and FE-SEM can satisfy the required standards for nanopore imaging, but they cannot provide 3D pore structures, which is crucial to obtaining the pore shapes and interconnected regions. FIB-SEM can obtain voxel dimensions of 10 nm and is also the only method to characterize the nanopores of the sample by 2D imaging and 3D spatial components, simultaneously. This method is time-consuming and expensive, but many studies have been performed on nanopore characterization using this method because of its accuracy. At present, FIB-SEM has become the mainstream tool for imaging nanopores (Keller et al. 2011; Chalmers et al. 2012a, b; Dewers et al. 2012; Bai et al. 2013). 3D images of shale are first obtained by the FIB-SEM device (Fig. 9) through shape correction, brightness correction, depth-of-field correction and phase distinguishing (Sun et al. 2016), which will make the results of the simulation more accurate and credible.

3D digital core of shale (Sun et al. 2016)

4.1.2 Knudsen flow

The Knudsen number (Kn) is the ratio of the gas molecular mean free path to the characteristic length. It is an important parameter to define the gas flow regime in shale. The gas transport though shale can be characterized by the value of Kn as continuum flow (Kn < 0.0001), slip flow (0.0001 < Kn < 0.1), transition flow (0.1 < Kn < 10), and free molecular flow (10 < Kn) (see Fig. 10) (Javadpour et al. 2007; Sondergeld et al. 2010; Wang et al. 2015). There is a deviation of the gas flow from the continuum hypothesis because of the influence of high Knudsen numbers in shale (Zhang et al. 2014), and other rarefied gas flow mechanisms begin to dominate the gas migration (Wang et al. 2015). Under this condition, Navier–Stokes (N–S) equations with no-slip boundaries will be invalid (Roy et al. 2003). Conventional numerical simulation methods based on the N–S equations also lack the foundation to correctly describe the transport characteristics of shale gas or other rarefied gases.

Relationship among Kn, pore diameter and mean pressure of reservoir at 400 K (Wang et al. 2015)

4.1.3 Adsorption and desorption

Shale contains much organic matter. Previous works suggest that CO₂ has stronger interactions with the organic constituents compared with shale gas because of the appreciable electrostatic moment of the CO₂ molecule (Middleton et al. 2015). Many experimental studies have also revealed that CO₂ can efficiently displace absorbed shale gas (Kang et al. 2011; Liu et al. 2013a, b; Heller et al. 2014; Luo et al. 2015). However, the displacement mechanisms of adsorbed CH₄ by CO₂ injection are far from being well understood, and quantitatively accurate assessments must be provided. For example, the popular Langmuir model neglects the interaction between CH₄ and CO₂ and does not consider nonidealities in the fluid phases (Middleton et al. 2015). Recently, Shi and Durucan (2003) developed a bidisperse pore-diffusion model for the displacement and desorption of CH₄ in coal by CO₂ injection. Wu et al. (2015) and Yuan et al. (2015) first employed molecular dynamics (MD) simulations to study the displacement of adsorbed CH₄ by CO₂ injection in carbon nanochannels, and they obtained the breakthrough time and front shape of the CO₂. However, to the best of our knowledge, there have been few related reports about the CO₂–CH₄ displacement structures in the slit pore. The adsorption and desorption of gas can induce slit pore changes and it provides feedback that impacts the gas flow. For example, the primary slit pore of shale narrows because of CO₂ adsorption-induced swelling (see Fig. 11, marked A) (Yin et al. 2016). Therefore, these issues cause great difficulties in the accurate simulation of the CO₂–CH₄ displacement process.

CO₂ adsorption-induced swelling for shale. a, c, e Untreated and b, d, f ScCO₂-treated (Yin et al. 2016)

4.2 LB-MD multiscale method

Molecular dynamics (MD) simulation is an effective method to study complex flows at the nanoscale (Li et al. 2010). In MD simulations, the fluid behaviour is described by the motion of the individual particles interacting with each other via intermolecular potentials (Koplik and Banavar 1995). As a result, it is only suitable to simulate fluid behaviours at the nanoscale owing to the limitation of computation cost (Zhao et al. 2016). LBM, as an intermediate method between continuum and atomistic simulations, is considered to have tremendous potential for simulating the fluid flow at the meso-scale, and detailed information on it has been presented in Sect. 2.2. However, the interaction force between the fluid and rock is assumed to be directly proportional to the density of the fluid in conventional LB models, which has no theory to support it and does not represent the real physical phenomena in certain situations (Liu et al. 2015a, b). To combine the advantages of the two described methods, coupling the LBM and MD is proposed to simulate fluid flows. MD describes the atomistic recognition of the fluid flow process near the walls of the porous media, while LBM is responsible for the rest of the simulation process and decreases the computational cost (Marsh et al. 2010). Combinations of the LB and MD methods can be divided into two categories. The first type is the synchronous coupling of the two methods performed by Ahlrichs and Dünweg (1998, 1999) and Fyta et al. (2006). In their studies, the simulation of the motion of DNA, polymers and other macromolecules in water was performed by the method. However, such a synchronous calculation method is very time-consuming for the simulation of fluid flow in porous media owing to the large computation cost of parameter passing between MD and LBM. The second type is the coupling of MD and LB in tandem to settle a problem, rather than simultaneous coupling (Marsh et al. 2010; Choi et al. 2014; Liu et al. 2015a, b; Pereira et al. 2016). MD provides the atomistic recognition of the flow physics near the walls, while LBM describes the rest of the simulation domain (including the Knudsen flow and absorption/dissolution-induced pore structure changes) and substantially reduces the computational cost.

To accurately simulate the CO₂-shale gas displacement process in the pore structure of shale, the 3D digital core of shale that is obtained by the FIB-SEM device is used directly for the LBM simulation with MD. Considering the computational cost, the coupling of MD and LB in tandem is employed to simulate the CO₂-shale gas displacement process (Fig. 12).

Schematic diagram of the workflow

First, the interaction forces between CO₂–CH₄ are obtained from the P–R EOS calculation (Liu et al. 2015a, b), and the shale matrix–CO₂/CH₄ interaction (including the displacement of adsorbed CH₄ by CO₂ and the interaction force of fluids-shale) is acquired from the MD simulation. According to the MD simulation, we can calculate the interfacial tension of CO₂/CH₄, the density of CO₂/CH₄ and the concentration of CO₂/CH₄. Then, we run the lattice Boltzmann simulation based on the 3D digital core of shale to obtain the real and whole CO₂–CH₄ displacement process.

CO₂-formation fluid displacement at the pore scale is a key issue in CO₂ sequestration and enhanced oil/gas recovery. The current situation and limitations of experiments and numerical simulations in studying the CO₂-formation fluid displacement at the pore scale are described in detail. From the review, we can conclude that there are many influencing factors for CO₂ displacement at the pore scale, e.g., displacement type, effective porosity, flow rate, Ca and M. The use of MRI or X-ray CT and the observation of the fluid displacement in fabricated micromodels are recognized as effective approaches to detect and observe the CO₂-formation fluid displacement processes. The LBM can precisely simulate the CO₂-formation fluid displacement processes in a 2D or 3D complex pore structure, but at the present, realistic rock pore geometries are not used directly to simulate the displacement process based on LBM. Various enhancements to overcome the limitations of experiments and numerical simulations in studying the CO₂-formation fluid displacement are suggested. Applications of CO₂ displacement have been highlighted for shale gas production. Therefore, the status of research and challenges in the application of CO₂ displacement to enhanced shale gas production are reviewed in the paper. Then, the coupling of MD and LB in tandem is proposed to simulate the CO₂-shale gas displacement process in the 3D images of shale obtained by an FIB-SEM device.