Abstract

Understanding fluids migration and leakage risk along the fault zone is necessary to guarantee the safety of CO₂ geological storage. The validity of Darcy’s law gets challenged in dealing with the flow in open fractures since the occurring of turbulence flow. In this study, we develop a 2D model with usage of T2Well, an integrated wellbore-reservoir simulator, to investigate the leakage problem along open fractures which are embedded in a fault zone from the deep injection reservoir to shallow aquifers. The results record a positive feedback of gas expansion and pressure response in fracture, which causes a quick downward propagation of highly gas saturated zone from the top of fracture and an easy gas breakthrough in the shallower aquifers. The decreasing of aperture size of fracture significantly enhances the leakage rates in fracture, but with less influences as aperture increases. In comparison, the Equivalent Porous Media models show a good approximation with the momentum model of large apertures but poor for the small one. Nevertheless, the differences are small in terms of final CO₂ distribution among various aquifers, suggesting that Darcy’s law may be still “effective” in solving flow problem along fractures in a constant injection system at a large time scale.

1. Introduction

Injection of CO₂ into deep saline aquifers for CO₂ geological storage (CCS) is considered as an effective technology to mitigate greenhouse gas emission. Leakage through existing fault into shallow groundwater aquifers, which occurs in many natural systems (e.g., [1–4]), is an important risk associated with CCS. Many numerical studies highlight this risk [4–9]. Pruess and Garcia [10] investigate CO₂ discharge along a fault zone by a simplified 1D model of constant boundary conditions, to find that fluids dynamics in the fault zone are mainly constrained by relative permeability and capillary function. Lu et al. [11] set up a 2D fault model and use time-varying boundary conditions to investigate the risk of CO₂ leakage through a fault. The results suggest that fault permeability is the most sensitive factor affecting both CO₂ and brine leakage rate. Some studies further investigate the impacts of spatial heterogeneity about hydraulic parameters distribution in the fault zone on CO₂ leakage dynamics. Vialle et al. [6] develop a 2D model of CO₂ leakage along a fractured cap rock, introduce heterogeneities in the initial permeability field of fractured damage zone, and find that the leakage of CO₂ into upper aquifer does not spread uniformly but most likely show one or several points’ source leaks. Rinaldi et al. [12] suggest that the presence of hydraulic heterogeneity in the fault would influence the pressure diffusion when investigating the effect of hydraulic properties variation within a fault zone embedded in a multilayer sedimentary system. Jeanne et al. [9] get similar conclusion when studying the impact of hydraulic-property variations on fluids flow along a fault zone by two damage-zone models. They find that CO₂ is trapped in the fault zone and obtain a high fluid overpressure if the heterogeneous model is used. In contrast, the homogeneous model predicts an easy migration of CO₂ with lower overpressure. By recognizing the fact that fracture is usually the major flow pathway in a fault, a few researchers recently introduce the discrete fracture network (DFN) model to study fluids flow through fault zone or fractured rocks [13–15]. In general, these models divide reservoir system into fracture network and matrix rocks and approximate a fracture as open space bounded by parallel plates. The equivalent permeability of a fracture is usually estimated using the cubic equation [16].

In the above studies, the fault is all assumed to be a kind of porous media so that the flow can be described by Darcy’s law even in the models with preferred flow path such as open fractures. However, a fracture is fundamentally an open space bounded by two walls. Under some conditions, the aperture could be larger enough so that the flow in the fracture becomes turbulence or the momentum term is no longer negligible compared with the friction term. In this situation, Darcy’s law is no longer proper to describe the flow dynamics and the full Navier-Stocks momentum equation would be needed. Furthermore, when free CO₂ phase evolves, the flow dynamics become more complicated because of complicated phase interference (e.g., gas-lifting phenomena) which may not be accounted by a Darcy law based model.

In this paper, we will present a numerical modeling study of CCS induced leakage from the deep injection aquifer to the shallow aquifer using T2Well [17], an integrated reservoir-wellbore simulator that solves the compressible two-phase Navier-Stocks momentum equations for flow in the subdomain of open space (e.g., well or fracture). The hypothetical leakage problem is formulated as a 2D problem (-plane) so that the flow in open fractures can be simplified as a 1D problem to facilitate the usage of T2Well. The focus is on the impacts of open fractures on the possible leakage flow dynamics and under what conditions the flow in the open fracture can be properly approximated as Darcy flow.

2. Method

2.1. Concept Model and Numerical Grid

The problem we are interested in is the possible leakage through a fault from a deep injection aquifer to the shallow ones due to the injection of CO₂ (Figure 1). We take a profile of the system with a thickness of 2 m in -direction as our model that can be simplified as a 2D leakage problem without losing generality (Figure 1). As shown in Figure 1, in the -plane, there are two sets of aquifers. The two shallow aquifers (Aquifers 3 and 4), separated by an aquitard in vertical direction, are open to constant boundary conditions on the left end but bounded by the fault at the right end. The two deep aquifers (Aquifers 1 and 2), also separated by an aquitard in vertical direction, are bounded by the fault at the left end. The fault zone is located at the distance of 4500 m from the left boundary and 500 m from the right boundary of the model. Aquifer 1 is the injection aquifer where a horizontal injection well is drilled at its right and bottom end for successive injection of CO₂ with a constant rate of 0.25 kg/s per unit meter for 5 years. As the pressure and CO₂ plume reach the fault, they may propagate further into the shallow aquifers (3 and 4). Aquifer 2 may serve as a buffer zone (or thief zone) in some degree depending on the pressure conditions at the nearby fault. The fault zone provides the critical pathway in this leakage problem. It is widely expected that the major flow pathway in a fault zone is the well-developed fracture network. Because flow modeling with detailed geometry of the fracture network is often not practical (if not possible), we approximate the major flow pathway as a virtual fracture. The virtual fracture can represent different fracture network (i.e., the major flow pathway) with its effective aperture. To do so, we first define a fixed fracture zone in the numerical grid (e.g., 5 cm width in -coordinate) and then use perimeter of the virtual fracture, Г, and the effective porosity of the fracture zone, , defined as the fraction of the area occupied by fracture and total cross-section area of the fracture zone. Figure 2 shows the sketch of the virtual fracture model. As shown in Figure 2, with being fully fractured in direction of the fracture zone, the porosity can be calculated as the ratio of fracture aperture, , over the width, , of facture zone. The fracture perimeter, , can be calculated as . For example, a fracture aperture of 1 cm width (Case  1 in Figure 2) will result in a porosity of 0.2 and perimeter of 402 cm. If the aperture decreases to be 1 mm (Case  2 in Figure 2), the corresponding porosity is 0.02 and perimeter 400.2 cm. Similarly, we can have a porosity of 1.0 and perimeter of 410 cm if the effective fracture aperture is 5 cm. In this way, the virtual fracture as represented as the fracture zone can be seen as a lumped flow pathway representing the entire fracture network in the fault.

Figure 1 

Concept model of CO₂ leakage along the fracture zone developed in the middle of a fault area which offsets and connects two sets of aquifers, the deep Aquifers 1 and 2 and shallow Aquifers 3 and 4. Open fractures of different apertures may grow inside the fracture zone and work as a fast pathway for CO₂ leakage. A horizontal well of 2 m length is drilled at the right and bottom side of Aquifer 1 for 5 years’ CO₂ injection.

Figure 2 

A cross-section (-plane) view of the fracture zone and relationships of the fracture aperture, , and perimeter, , and the porosity, , of the fracture zone in different cases of various aperture sizes. Note that the cross-section of the fracture zone remains unchanged in our numerical grid (2 m in length, , and 5 cm in width, ) but the size of the virtual fracture (the blue colored area) can change and is defined by the perimeter of fracture space and the effective porosity of fracture zone and the ratio of fracture space over the total cross-sectional area of the fracture zone. In picture, the block-shaped area denotes the solid impermeable rocks in the fracture zone. In Case  1, a fracture of 1 cm wide aperture would have a perimeter of 402 cm (2 × 200 + 2 × 1 cm) and a porosity of 0.2 (1 × 200 cm²/5 × 200 cm²). The small aperture case ( mm) shown in Case  2 results in a porosity of 0.02 and perimeter of 400.2 cm. In particular, in Case  3, where a fracture fully develops among the fracture zone, its aperture finally equals the width of fracture zone of 5 cm, and the porosity becomes 1.0.

Figure 3 represents the mesh we construct based on the concept model. The dimension extends from −4500 m to 500 m in -axis and −1500 m to −200 m in -axis. The fault zone located at the zero in -axis is composed of a 0.05 m wide fracture zone and several surrounded damage zones of width closing to 1 m. The horizontal resolution of the grid varies from 0.05 m near the fracture to the maximum of 50 m and 500 m at the right and left boundary, respectively. Table 1 gives all hydraulic parameters of the porous medias used in this study while the information about fracture zone is summarized in Tables 2 and 3. The rock properties of porous media derive from the wellbores drilled at the Yaojia formation in SanZhao depression of Songliao Basin, China. Commonly distributed sandstone strata with high porosity and permeability and the coupled upper shale of poor hydraulic property in Yaojia formation make the reservoir suitable for oil and CO₂ storage. In this study, we select a large porosity and corresponding permeability for the aquifers and small ones for aquitard, clays, and cap rocks, to guarantee the most of injected CO₂ can be sequestrated among aquifers instead of other media.

(a)

(b)

(a)
(b)

Figure 3 

Mesh discretization based on the concept model shown in Figure 1. Picture (a) shows a full view of the grid in which a CO₂ injection well (the cyan star) is located at the right bottom side of Aquifer 1. Picture (b) is the zoom view of the fault zone, comprising a 0.05 m wide fracture zone in the middle and (highlighted by the cyan solid line) several surrounded damage zones of width closing to 1 m. The observation points B1, B2, and B3 are located at = −1290, −975, and −410 m, respectively ( = 0.0 m).

Before CO₂ injection, the whole system is saturated with water and under hydrostatics pressure. Surface temperature is 20°C and the temperature gradient is set to be 3°C/100 m. As shown in Figure 3, three points B1, B2, and B3 are located at the bottom ( m), middle ( m), and top ( m) of the fracture to monitor the pressure responses there associated with CO₂ injection. In addition, the fluids status and flow dynamics at injection well and the interfaces between the fracture and each aquifer are also monitored.

2.2. The Numerical Simulator

We use T2Well/ECO₂N [17] to simulate the fluid and heat flow associated with the leakage problems in this study. T2Well is an integrated reservoir-wellbore simulator. The major governing equations used in T2Well are listed in Table 4. As shown in Table 2, for the subdomain of porous media (e.g., formations), the multiple phase version of Darcy’s law is used to obtain the phase velocity whereas the phase velocities are calculated as a function of the mixture velocity and the drift velocity in the subdomain of wellbore. By applying the empirical Drift Flux Model (DFM), the two-phase momentum equations can be simplified as a single momentum equation in terms of the mixture velocity:where is the cross-sectional area of the wellbore and the term is caused by slip between the two phases. The terms , , and are the mixture density (), the mixture velocity (), and the profile-adjusted average density of the mixture (). is the perimeter of the cross-section, and is apparent friction coefficient, a function of wall roughness and flow regime.

The details of the approach used in T2Well for modeling flow in wellbore can be found in the manual or other related papers [17–19] and would not be duplicated here. However, we have slightly modified the code to adapt the simulator to the case of fracture flow because the open space of a fracture is not of a circular shape as in the case of a wellbore. The momentum equation (1) uses the general terms, the perimeter Г, and the cross-section area , which should have no problem to apply to any geometric shape of open space including the rectangle shape of a fracture. The major difference is in the friction coefficient, , which is calculated as a function of local Reynold number, the wall roughness, and the effective diameter of the pipe in T2Well [18, 19]:where is the roughness of the wellbore and the Reynold number Re is defined aswhere is the mixture viscosity. In both (2) and (3), the effective diameter can be seen as the critical length for a pipe. Therefore, the modification we made is to let T2Well use the aperture in the place of the effective diameter in calculation of the Reynold number and the friction coefficient for the subdomain of open fractures.

ECO₂N is a TOUGH2 EOS module for modeling thermal physical behaviors of the H₂O-CO₂-NaCl system (Pruess, 2005) and has been widely used in modeling flow problems associated with CCS. In this study, we use an updated version of ECO₂N (i.e., Version 2.0, Pan et al., 2017) which provides more accurate calculations of the thermal physical properties of the CO₂-rich phase and can be applied to higher temperature range than the previous version (i.e., Version 1.0). Following the convention of ECO₂N, we will call the CO₂-rich phase gas phase thereafter although it could actually be gaseous, liquid, or supercritical CO₂ with dissolved water depending on local pressure and temperature.

3. Results and Discussions

3.1. CO₂ Migration Dynamics in Case  1 (Aperture = 1 cm)

Injection of CO₂ increases the pressure at the injection point first (Figure 4). The pressure propagates quickly from injection point to the left end of Aquifer 1 which causes water to start flowing out from Aquifer 1 to the fault at about 0.1 days (Figure 5). The pressure responses at different depth in the open fracture (i.e., at Points B1 through B3) are almost identical indicating that the pressure propagation through the open fracture is much faster than through the porous formation (Figure 4). However, the responses in terms of gas phase saturation are very different. The gas saturation at shallow depth in the fracture can be higher than that at the injection well (Figures 4(a) and 4(d)) while the gas saturation at B1 (the deepest point in the fracture) is still below 0.012 at end of simulation (Figure 4(b)).

(a)

(b)

(c)

(d)

(a)
(b)
(c)
(d)

Figure 4 

Pressure buildup and gas saturation changes at injection well and monitor points during CO₂ injection in Case  1.

Figure 5 

Mass flow rate evolution of the gas/liquid phase and total CO₂ at the outlet of Aquifer 1 in Case  1.

After about 4 days, CO₂ starts to flow into the fault as dissolved CO₂ in water (Figure 5). This period lasts for 254 days before breakthrough of gas phase. During this period, the CO₂-dissolved water flows up along the fracture and quickly reaches the top of fracture. As a result, the mass fraction of the dissolved CO₂ increases with time but its contour line is close to vertical above Aquifer 1 (Figure 6(a)).

(a)

(b)

(c)

(d)

(a)
(b)
(c)
(d)

Figure 6 

Distributions of the dissolved CO₂ mass fraction , gas saturation , and gas phase flow rate along the fracture profile in Case  1 during the few days around the breakthrough of CO₂-rich phase at 258th day. Note that the fraction of dissolved CO₂ constantly becomes saturated if gas phase CO₂ evolves, working as a flag to trace the occurrence of two phases in fracture.

Gas phase occurs first at bottom of the open fracture at day 257.9 although the gas saturation is very low (Figure 6(b)). The bubbling region expands to the depth of 700 m almost immediately although in much lower overall gas saturation at day 258.2. After that, the gas saturation in the region gradually increases until a more significant increase occurs in the neighborhood of the 730 m depth where CO₂ transforms from supercritical to gaseous phase, which results in a large expansion. This period ends at day 258.5 as the bubble region suddenly reaches the top of the open fracture and, as a result, the overall gas saturation below 700 m depth decreases. The occurring of gas phase at top of the open fracture, where the pressure is the lowest in the fracture, starts positive feedback that results in a quick expansion of high gas saturation region downward (Figures 6(b) and 6(d)). This positive feedback can be described as follows. With increase of gas phase saturation at top of the fracture, the overall density of fluid mixture decreases significantly because the density of gas phase is much smaller than the aqueous phase. As a result, the gravity applied on the fluid below reduces greatly. This will in turn greatly reduce the pressure there which will cause more dissolved CO₂ to escape from the aqueous phase and the gas phase to occupy more volume. This process can go on until the difference in density between the gas phase and liquid phase is not large enough to make large change in gravity force due to change in gas saturation.

This downward propagation of high gas phase saturation region in the fracture has great impacts on the sequence of CO₂ breakthrough into each aquifer (Figure 7). The gas phase flow into aquifer occurs first at the shallowest aquifer (Aquifer 4) and then Aquifer 3 (deeper than Aquifer 4) about 11 days later (except a short lived inflow at day 258.5). For much deeper aquifer (Aquifer 2), there is no significant gas flow although it is the first aquifer in the leakage pathway after Aquifer 1. This is mainly because the gas saturation at that depth is too low (Figure 6(d)).

(a)

(b)

(c)

(d)

(a)
(b)
(c)
(d)

Figure 7 

Phases mass flow rate at the interfaces of fracture and various aquifers in Case  1 before and after gaseous CO₂ breakthrough.

The change in fracture pressure associated with the quick expansion of the high gas saturation region in the fracture shows great impacts on the flow dynamics between the fracture and the aquifers since the fracture is opened to multiple aquifers at various depths. Although the gas flow from Aquifer 1 to the fracture increases with time and the liquid flow decreases in general and then tends to stabilize at certain level, such increase or decrease occurs at step-like way associated with some oscillations (Figure 7(a)). Such oscillation is more profound at shallower depths (Figures 7(c) and 7(d)) and can even be found in the liquid flow rate between the fracture and Aquifer 2 (Figure 7(b)).

The profile of Reynold number shows that the flow in fracture is still mostly laminar flow with the maximum Reynold number < 2400 in the fracture (Figure 8(a)). The higher Reynold number region almost overlaps the gas phase dominant region while the Reynold number decreases with depth mainly because the velocity decreases with depth due to fluid leaking to Aquifers 3 and 4 on its way to top (Figures 8(a) and 8(c)). However, there are small regions where the momentum term is important because the ratio of the momentum over friction term is closing to or excesses 1.0.

(a)

(b)

(c)

(a)
(b)
(c)

Figure 8 

Profiles of Reynold number and the ratio of momentum over friction gradient at top zone of the fracture in Case  1 before and after gas phase evolves.

As shown in Table 5, after one year’s injection, 77% of the injected CO₂ is still in Aquifer 1 and about 20% leaks into shallow aquifers, among which Aquifer 4 shows much potential to steal CO₂ than Aquifer 3. However, at end of 5 years’ injection, only 20.2% of the injected CO₂ stays in Aquifer 1 while 42.5% leaks into Aquifer 3 and 35.3% into Aquifer 4 (Table 6). The time dependent characteristics of CO₂ mass distribution in each aquifer reflects the impacts of fluids dynamics at different time scale. Finally, there is no gas phase CO₂ sequestrated in Aquifer 2 with the increase of injection duration, although the dissolved CO₂ slightly increases.

3.2. Impact of Fracture Aperture

In a real fault, the fracture network development could vary widely depending on local stress field and other factors. We did a series of simulations of the same leakage process using the fracture aperture as the parameter characteristic of the given open fracture system developed in a fault. Two additional cases with smaller (Case  2, 1 mm) or larger (Case  3, 5 cm) apertures, shown in Figure 2 and Table 2, are simulated to investigate how the aperture affects the CO₂ migration along fracture, in terms of pressure responses, gas saturation, fluids dynamics, and the final CO₂ distributions in aquifers. The smaller aperture case represents a scenario with a poor fracture network developed in the fault whereas the larger aperture case represents the other end. As shown in Figure 9(a), the injection pressure quickly increases with time in all cases until the gas phase reaches the fault. After that point, the injection pressure increases much slowly because the leaking fluid accesses the lower resistance pathway (i.e., the open fracture). Different from the other cases, the injection pressure in the small aperture case (Case  2) significantly drops after that point which is companioned with the similar pressure drop at bottom of the fracture and very fast pressure increase at top of the fracture (Figures 9(b) and 9(c)). In view of overall pressure response in the fracture profile shown in Figures 9(e) and 9(f), Case  2 is also different from other cases that the pressure greatly increases at top and meanwhile decreases at bottom of the fracture. This is most likely because the gas saturation in fracture in Case  2 is much higher than other cases (Figures 10 and 6(d)) so that the resistance to gas flow is much smaller in Case  2 than in the other cases (Figure 11).

(a)

(b)

(c)

(d)

(e)

(f)

(a)
(b)
(c)
(d)
(e)
(f)

Figure 9 

Pressure responses of CO₂ breakthrough at injection well and the fracture.

(a)

(b)

(a)
(b)

Figure 10 

Gas saturation changes in the fracture in Case  2 and Case  3 associated with CO₂ breakthrough.

(a)

(b)

(c)

(a)
(b)
(c)

Figure 11 

Distributions of gas velocity in the fracture profile for Cases  1 to 3 as gaseous CO₂ evolves.

On the other hand, the fracture aperture shows great impacts on the flow dynamics in the fracture after two-phase flow has been established. The Reynold number increases as the aperture decreases from 1 cm (Case  1) to 1 mm (Case  2) because of increase in velocity (Figures 12(a) and 8(a)). However, the increase is quite small compared with Case  1 so that the maximum Reynold number is still below 2400 to prevail a laminar flow in the fracture. What is different is that the ratio of momentum term over friction becomes quite small as aperture decreases from 1 cm (Case  1) to 1 mm (Case  2), suggesting that the friction of side wall of fracture dominates fluids flow in the fracture other than the momentum term. In contrast, the Reynold number significantly increases over 2400 accompanied with momentum/friction ratio greater than 1.0 when aperture increases from 1 cm (Case  1) to 5 cm (Case  3). At these conditions, the fluids flow is of high turbulence and the momentum term becomes important in controlling the fluids flow rate in the fracture.

(a)

(b)

(c)

(d)

(a)
(b)
(c)
(d)

Figure 12 

Changes of Reynold number and momentum/friction ratio in the fracture of Case  2 and Case  3 when CO₂ breakthrough occurs.

The larger the fracture aperture is, the earlier the gas phase arrives to the fracture from the injection well (Figure 13(a)) mainly because larger aperture means less resistance to fluid flow. The difference is quite small when the aperture increases from 1 cm (Case  1) to 5 cm (Case  3) though. In contrast, when the aperture reduces from 1 cm (Case  1) to 1 mm (Case  2), occurring of the gas phase flow from Aquifer 1 to the fracture delays about 5 days but the mass flow rate quickly exceeds those in the case with larger aperture after gas phase reaches top of the fracture. This is consistent with pressure response pattern at fracture bottom (Figure 9(b)) that the pressure is much higher in Case  2 than in the other cases before gas phase breakthrough and the pressure in Case  2 after that point quickly drops to below those in other cases. As a result, the liquid phase flow rate from Aquifer 1 associated with CO₂ breakthrough drops slowly to be higher values than other cases in Case  2 (Figure 13(b)). The similar trend in terms of gas flow rate into shallow aquifers (Aquifer 3 and 4) can be found for aperture reduced to 1 mm, that is, Case  2 (Figures 13(c) and 13(d)). However, the occurring of gas flow into shallow aquifers delays as the aperture increases from 1 cm (Case  1) to 5 cm (Case  3) although the gas flow rate quickly approaches similar values in these two cases.

(a)

(b)

(c)

(d)

(a)
(b)
(c)
(d)

Figure 13 

Fluids dynamics at the outlet/inlet of each aquifer in Cases  1–3 at early time (up to 300 days) after CO₂ breakthrough.

In general, the liquid phase mass flow rate decreases with time as the gas phase mass flow rate increases in all cases (Figure 14). The exemption is Case  2 (1 mm) where much higher liquid flow rate into Aquifer 4 (including two peaks that exceed the rate during single liquid phase period) is found. This could be because the liquid is being lifted by the higher gas velocity in this smaller aperture case. The gas phase flow into Aquifer 2 is not significant in all cases probably because of very low gas saturation in the fracture at that depth. However, the liquid flow rate into Aquifer 2 shows different characteristics as aperture decreases from 1 cm to 1 mm (Figure 15). In general, the liquid flow rate is higher in the small aperture case (Case  2) before gas CO₂ enters the fracture especially at the early injection stage (e.g., 50 days), which is also supposed as the reason for the higher amount of dissolved CO₂ stored in Aquifer 2 at end of 5 years’ injection (Figure 16). With gas CO₂ breakthrough, the liquid reversely flows out from Aquifer 2 in Case  2 while other cases fail to observe this phenomenon. This is most likely because of the sudden pressure loss occurring at the fracture nearby Aquifer 2 as a result of the quick expansion of gas bubble zone at top of the fracture.

Figure 14 

Fluids dynamics in Cases  1 to 3 at the interface of Aquifer 4 and the fracture with CO₂ breakthrough.

Figure 15 

Evolution of liquid phase mass flow rate into Aquifer 2 from the beginning of CO₂ injection until the end of 300th days.

(a)

(b)

(c)

(d)

(a)
(b)
(c)
(d)

Figure 16 

CO₂ mass distribution into each aquifer at the end of 5 year’s injection in percentage.

Despite many differences in flow dynamics caused by the fracture aperture during the first year of injection, especially during phase transition period in the fracture, the final distributions of the injected CO₂ are very similar for all cases (Figure 16). This is because the gas flow rates into each aquifer approach the same values after first year (Figure 17), indicating that the fracture aperture only affects the flow dynamics at early time, especially during the phase transition period in the fracture, but it has insignificant effects on the quasi-steady state leakage corresponding to a fixed injection rate adapted in this study.

(a)

(b)

(c)

(d)

(a)
(b)
(c)
(d)

Figure 17 

Fluids dynamics at the outlet/inlet of each aquifer in Cases  1–3 for whole time scale of CO₂ injection.

3.3. Equivalent Porous Media (EPM) Approximation

In many modeling studies of leakage through open fractures in a fault, the EPM approximation is often used that the Darcy law is assumed to be valid for describing the flow dynamics in open fractures. We have made a few EPM models of the same problem with the effective permeability of the open fractures, being estimated as a function of the fracture aperture, , suggested by Zimmerman and Bodvarsson [20] and Witherspoon et al., 1980:The effective permeability of fracture zone, , can then be calculated as follows:where is the “porosity” of fracture zone, defined as the fraction of void space in fracture zone and is the permeability of matrix rocks surrounding the virtual fracture. In this study, a value of zero is preferable for since the virtual fracture is the only pathway for CO₂ leakage along the fracture zone. The calculated effective permeability of the fracture zone for Cases  4, 5, and 6 is shown in Table 3. Cases  4, 5, and 6 are corresponding to Cases  1, 2, and 3, respectively.

Figure 18 shows the differences of pressure response at the injection well and monitor points in the fracture for two pairs’ cases (Case  1 versus Case  4 and Case  2 versus Case  5) of the same aperture but different models. The pressure in large aperture cases ( cm in Case  1 versus Case  4) is similar although the EPM model records a bit higher pressure at injection well and the bottom of fracture after CO₂ breakthrough. With aperture increases from 1 cm to 5 cm, the pressure responses in Case  3 and Case  6 are almost identical (though they are not shown here to avoid too much crowd). However, when it comes to the small aperture cases ( mm in Case  2 versus Case  5), big differences are observed that the injection pressure drops slowly in EPM model but fast in momentum model after CO₂ breakthrough; the pressure of EPM model before CO₂ breakthrough at the bottom of fracture (B1 point) is much lower and drops less in response with the CO₂ breakthrough as compared with the momentum model; what is more, the pressure in EPM model gets less increase when gas CO₂ reaches the top of fracture (B3 point).

(a)

(b)

(c)

(a)
(b)
(c)

Figure 18 

Models comparison (Case  1  versus Case  4, Case  2 versus Case  5) in terms of pressure response at injection well and the bottom and top of fracture.

These differences of pressure response further determine the various characteristics in terms of fluids flow dynamics at each aquifer/fracture interface (Figure 19). Overall, the larger aperture cases ( cm in Case  1 versus Case  5) show good approximations between two models while showing bad ones for the small aperture cases ( mm in Case  2 versus Case  4). For instance, Case  2 records a much higher gas and liquid phase flow rate from Aquifer 1 and gas phase leak rate into Aquifer 3 and Aquifer 4 when gas phase arrives at the top of fracture, in contrast with the small values predicted by EPM model (Case  4).

(a)

(b)

(c)

(d)

(a)
(b)
(c)
(d)

Figure 19 

Differences of fluids dynamics at the outlet/inlet of each aquifer between momentum model and corresponding EPM model of various apertures.

These seemingly strange results that the flow in larger fracture can be modeled by EPM better than in smaller fracture are for reasons. In general, larger fracture has less resistance to flow so as to be easy to develop high velocity flow regime which often causes turbulence flow and complex two-phase flow phenomena which may not be modeled properly by Darcy’s law-based EPM model. However, in this particular study, the fracture is only the middle pathway between porous aquifers and the mass flow rate is limited by the fixed injection rate. Under this condition, the large aperture will result in smaller velocity because of larger cross-section area or volume of the fracture (Figure 11). As a result, the gas velocity in the larger aperture cases (e.g., 1 cm and 5 cm) is too small to enhance the liquid flow in the fracture for the given injection rate adapted in this study. However, when the aperture is reduced to 1 mm (Case  2), the gas velocity is high enough (Figure 11) to enhance the liquid flow in the fracture (Figure 19). This so-called gas-lifting effect cannot be captured by a Darcy’s law-based EPM model even though the overall resistance of the fracture to flow can be described by the effective permeability. Therefore, whether or not the EMP model is a good approximation not only depends on the geometry of the fractures but also depends on the flow conditions in a particular system. Nevertheless, from the view of CO₂ distribution among various aquifers at the end of 5 years’ injection, the differences become quite small with aperture changes, suggesting that Darcy law may be still “effective” in solving fluids flow along fractures at large time scale, especially when the fractures are only an intermediate pathway between two porous media formations like the one studied here.

4. Conclusions

Fault area is one of the most concerns in conducting CCS projects for the safety of storage. Complex network of fractures developed in a fault area as a result of tectonic or anthropogenic activities may work as the fast pathway for CO₂ migration to favor a channel flow especially when two-phase flow occurs. In this study, we set up a 2D model of CO₂ leakage through a fault zone (simulated as a virtual vertical fracture) which connects two sets of aquifers at different depth. The momentum equation is applied for the calculation of fluids velocity in the fracture, with the DFM model to handle phase interactions in possible two-phase flow. Our motivation is to figure out what the differences are about channel flow and Darcy flow in open fractures and under what conditions the flow in the open fracture can be properly approximated as Darcy flow.

The results indicate that injection of CO₂ quickly increases the pressure in reservoir to facilitate the migration of water and gas into the fracture. When gas enters the fracture, it firstly evolves at the bottom at very low saturation and then quickly flows upward. Once it reaches the top, the bubble expands quickly because of low pressure which increases the local gas saturation and reduces the average density of the gas-liquid mixture. This will lower the pressure below and causes the increases of gas saturation there due to escape of the dissolved gas and expansion of the gas phase. Such positive feedback will propagate downward until the density difference between the gas phase and the aqueous phase becomes small or the local mass fraction of gas component becomes very small. As a result, the occurring of gas phase flow into Aquifer 4 is earlier than Aquifer 3 although Aquifer 3 is below Aquifer 4. With aperture decreases from 1 cm to 1 mm, the friction force increases and, as a result, the momentum/friction ratio sharply decreases, although the Reynold number gets less influenced and remains below 2400 to prevail the laminar flow in the fracture. In contrast, when aperture increases from 1 cm to 5 cm, the Reynold number and momentum/friction ratio significantly increase to be the maximum of 15000 and 10, respectively, suggesting the occurrence of turbulent flow. In case of small aperture (1 mm), the gas phase velocity is high enough to enhance the liquid flow in the fracture by gas-lifting effects. As a result, the fracture dries quickly and the pressure drops quickly after two-phase flow developed in the fracture in this small aperture case. The differences of pressure responses further lead to diverse gas saturation profiles along fracture as well as the fluids dynamics of each aquifer. At the end of 5 years’ injection, more dissolved CO₂ is stored in Aquifer 2 of Case  2.

With the relative small injection rate adapted in this study, the Equivalent Porous Media (EPM) model was found to be a good approximation of open fracture flow in cases of larger apertures but becomes poor if the aperture is reduced to 1 mm. The major reason is that, in the latter case, the gas velocity is high enough to enhance the liquid flow (i.e., the gas lifting) in the fracture which cannot be simulated by the EPM. Therefore, whether or not the EMP model is a good approximation not only depends on the geometry of the fractures but also depends on the flow conditions in a particular system.

In terms of CO₂ distribution among various aquifers at the end of 5 years’ injection, the differences due to different aperture is quite small except for the dissolved CO₂ in Aquifer 2. This is because the open fracture is only the intermediate pathway sandwiched by porous aquifers. Even though the early time flow dynamics, especially during the phase transition period, could be very different, the long term quasi-steady state flow would be quite similar for the given constant injection rate because the leakage is finally controlled by the high resistant media (i.e., the aquifers).

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The authors gratefully acknowledge funding by the Project of National Sciences Foundation of China (no. 41172217).