On Fluid and Thermal Dynamics in a Heterogeneous CO2 Plume Geothermal Reservoir

CO2 is now considered as a novel heat transmission fluid to extract geothermal energy. It can achieve both the energy exploitation and CO2 geological sequestration. The migration pathway and the process of fluid flow within the reservoirs affect significantly a CO2 plume geothermal (CPG) system. In this study, we built three-dimensional wellbore-reservoir coupledmodels using geological and geothermal conditions of Qingshankou Formation in Songliao Basin, China. The performance of the CPG system is evaluated in terms of the temperature, CO2 plume distribution, flow rate of production fluid, heat extraction rate, and storage of CO2. For obtaining a deeper understanding of CO2-geothermal system under realistic conditions, heterogeneity of reservoir’s hydrological properties (in terms of permeability and porosity) is taken into account. Due to the fortissimo mobility of CO2, as long as a highly permeable zone exists between the twowells, it ismore likely to flow through the highly permeable zone to reach the productionwell, even though the flowpath is longer.Thepreferential flow shortens circulation time and reduces heat-exchange area, probably leading to early thermal breakthrough, which makes the production fluid temperature decrease rapidly. The analyses of flow dynamics of CO2-water fluid and heat may be useful for future design of a CO2-based geothermal development system.


Introduction
The enhanced geothermal system (EGS) is defined as an engineered reservoir that has been created to extract economical amounts of heat from geothermal resources of low permeability and/or porosity [1].As part of an effort to reduce atmospheric emissions of carbon dioxide (CO 2 ), a novel concept of operating the EGS using CO 2 instead of water as the working fluid (CO 2 -EGS) and achieving simultaneous geologic sequestration of CO 2 has been proposed and evaluated [2,3].
In recent years, a similar concept, the so-called CO 2 plume geothermal (CPG) system, has been proposed [4].The CPG system utilizes existing, naturally porous, high permeability geologic formations (reservoirs) for geothermal energy recovery.The major benefit of the CPG system over the EGS is that the CPG system does not require hydrofracturing, which helps increase fracture permeability but may induce seismicity.The EGS has encountered considerable unfavorable conditions and sociopolitical issues (resistances).Consequently, the CPG system that can use the CO 2 sequestration site to recover geothermal energy may be practical.
The advantages of carbon dioxide as a working fluid compared to water in geothermal energy recovery include (1) large expansivity and compressibility, which result in great density difference between the injector and the producer and reduce the power consumption for circulation on account of buoyancy force, (2) low viscosity coefficient, which leads to the larger flow rate under fixed differential pressure, (3) low salt solubility, and (4) low chemical reactivity; CO 2 does not mostly tend to react with rocks.Pruess [3] has built a homogeneous five-spot "fully developed" fractures reservoir model and indicated that, for a given pressure difference between injection and production wells, CO 2 would generate 50% larger net heat extraction rate and four times larger mass flow rate compared to water.Atrens et al. [5] drew a conclusion without consideration of frictional pressure that CO 2 thermosiphon could produce similar amount of electric power with the same quantity of heat extraction compared to water but with simpler surface equipment.Xu et al. [6] performed wellbore-reservoir coupled five-spot model and pointed out that the specific enthalpy change of CO 2 in the wellbore is so small relative to its intrinsic specific enthalpy that the wellbore flow can be regarded as isenthalpic process, and the temperature of the fluids in production well will decrease rapidly with pressure.Compared to water, the advantage of CO 2 as a working fluid is more noticeable in low permeability and low temperature reservoirs.
Actually, the investigations mentioned above were carried out under an assumption that the target reservoirs are homogeneous.However, the permeability and porosity heterogeneity of a reservoir significantly affect the CO 2 migration and production.Garapati et al. [7] analyzed the effects of multilayered reservoirs on a CPG system; they found that the produced CO 2 mass fraction is dominated by high permeability layers and their position within reservoirs.In addition, the heat extraction rate comes down as the permeability of the bottom reservoirs decreases.Yang et al. [8] considered the effects of permeability and porosity heterogeneity on the liquid invasion in tight gas reservoir and drew a conclusion that both the liquid invasion depth and invasion rate increased as the heterogeneity coefficient of permeability decreased.Tian et al. [9] established a twodimensional model with the consideration of heterogeneity in hydrological parameters of a caprock and indicated that the hydrological heterogeneity significantly affects the containment of intruded CO 2 within a caprock.Bu et al. [10] constructed a 2D numerical model to study the influence of high-porosity and high permeability faults within a reservoir and found that the high-porosity and high permeability faults influence CO 2 migration and spatial distribution and then increase the CO 2 storage capacity of reservoir.
In this paper, we built a three-dimensional geological model with spatially varying permeability and porosity, based on Qingshankou Formation in Songliao Basin, China, to evaluate the effects of hydrological heterogeneity in terms of permeability and porosity on the performance of a CPG system.The temperature change, temporal and spatial distribution of CO 2 , rate of production flow, heat extraction rate of system, and storage capacity of CO 2 were chosen as the metrics to evaluate the performance.

Problem Setup
A great deal of detailed information is required to assess the feasibility of injected CO 2 as a heat transmission fluid at any specific site and to develop engineering designs for CO 2based geothermal systems.Before moving into site-specific investigations, general features and issues should be explored.This can be done by investigating deep brine systems that abstract site-specific features and thereby attempt to represent characteristics that are common to many such systems.Here, geological characteristics and thermophysical conditions are mainly extracted from the central depression of Songliao Basin, Northeastern China.The basin has a pretty high geothermal gradient and heat flow among sedimentary basins in China and can meet the temperature requirement for geothermal development.

Conceptual Model and Boundary Condition Settings.
To investigate processes of fluid migration and heat exchange within both wellbores and geologic formation, a threedimensional conceptual model with a size of 10000 m × 10000 m × 100 m, including two wellbores and a reservoir formation (Figure 1), to represent a 100 m thick sandstone layer extending 10 km horizontally in Qingshankou Formation in Songliao Basin, was used.The distance between the injection and production wellbores is 600 m, and both wellbores vertically perforate the whole thickness of reservoir.The length of the model set to be 10000 m laterally is aimed at avoiding the influence of lateral boundaries.Moreover, for the lateral boundaries, Dirichlet-type condition (here, constant pressure and temperature) was imposed.As for the top and bottom boundaries, no-flow conditions were assigned (see the following section for more details).The wellbores with the surrounding formation are also set to be no-flow conditions but the heat exchange is allowed.The mesh of the modeling domain is generated using TOUGHVISUAL [11], which is able to create regular and irregular grids for TOUGH family codes [12].During simulation, the accuracy and computational demand should be considered.Therefore, the mesh for the area around wellbores is refined, as shown in Figure 1.

Geological and Thermophysical
Parameters.The top of Qingshankou Formation in Songliao Basin is at 2510 m depth and the thickness of this formation is about 250∼550 m, including about 100 m thick sandstone layer with porosity of 10%∼30% and permeability of 3 to 300 millidarcys which is the main oil and gas reservoir in this formation and with overlying and underlying mudstone as the caprock and bedrock [13].The average geothermal gradient of the basin is 38.7 ∘ C/km [14].In our model, the depth of the sandstone reservoir is 3000 m with the thickness of 100 m.The initial temperature of the reservoir was set to be 120 ∘ C according to the average geothermal gradient and the local mean annual surface temperature (about 4 ∘ C).The initial pressure of the model was obtained by hydrostatic equilibrium, and the reservoir was initially saturated with resident water (i.e., CO 2 saturation Sg = 0).The details of reservoir geological and thermophysical conditions used in our model are listed in Table 1, and the specification of physical parameters of the wellbores is summarized in Table 2.

Heterogeneity Implementation.
From a realistic view, natural aquifer is intrinsically heterogeneous.The fluids are more likely to flow through the highly permeable zone, especially for CO 2 , because of its strong mobility.It may lead to an early breakthrough.Among the properties of reservoir, permeability and porosity have the largest influence on flow field and heat extraction.Therefore, we just take them as the main variables to study heterogeneity effects.
As shown in Figure 1, the reservoir extends 10 km laterally, including a one square kilometer subdomain in the center.Here we assume that the subdomain is a geothermal reservoir with heterogeneities in permeability () and porosity ().Neuzil [15] observed that there exists a log-linear relationship between permeability and porosity in argillaceous sediments.Tian et al. [9] drew a similar conclusion by regression analysis.In addition, Nelson [16] suggested a linear relationship between log() and  in quartz sandstone.In this study, we assumed that the above discussed relationship between permeability and porosity is suitable for the Qingshankou Formation, and the relationship can be expressed as where  and  are data-fitted parameters,  is the permeability in m 2 , and  is the porosity.Based on the data in Table 1, we figured out the value of  and  in our model as 10 and −15.52, respectively.We assumed that the permeability of reservoir is subject to a lognormal distribution, and thereby the porosity obeys a normal distribution inferred from (1).More realistically, the permeability obeys spatial random distribution in the reservoir rather than simply random distribution in probability and statistics.So we introduced the variation function, a concept derived from geostatistics, to depict the spatial distributions of permeability and porosity in the reservoir.In this study, a variance of 0.3 and correlation length of 300 m were chosen [9].The simulator T2well provides a flexible function inherited from TOUGH2 V2 for permeability modification for each individual grid expressed as where   is the absolute permeability of grid , as specified in data block ROCKS and   is the permeability modification coefficient which can be defined internally or externally and may be provided as part of the geometry data in block ELEME and then employed to multiply the absolute permeability   for each grid in the subdomain.When the permeability modification is accomplished, the porosity is yielded according to (1) and stored in block INCON.More details about permeability modification coefficient generation and heterogeneity realization can be found in previous research [9].The permeability and porosity range is shown in Table 3. Infinite types of heterogeneous permeability and porosity fields can be generated randomly.As shown in Figure 2, six cases with heterogeneity in permeability and porosity were selected to represent the conditions: (i) high permeability belt connecting two wellbores (hh1, hh2); (ii) highly permeable zones expanding away from the production wellbore (hp1, hp2); and (iii) low permeability zone between the two wellbores (hl1, hl2).They are not expected to match with any existing geologic formation; we aim to make some theoretical analysis and quantitatively evaluation of the effects of  heterogeneities on a CPG system and elucidate the affecting mechanism.

Modeling Scenario.
The CO 2 migrates simultaneously as the CO 2 is injected into the reservoir through injection wellbore, which pushes the existing resident water to move.Therefore, the production fluid is expected to be water at the early age.So the injection and production pressures at the wellhead are set to be 12 MPa and 0.1 MPa (atmospheric condition) to make CO 2 get breakthrough as soon as possible.
After a period of propagation, CO 2 passes through the reservoir and can be observed from the production wellhead, and then the saturation of CO 2 (Sg) in production fluids increases rapidly, and the temperature of carbon dioxide decreases sharply with the pressure decline in the production wellbore, which is called the Joule-Thomson effect.To get a higher production fluid temperature and control the flow rate for a stable operation of the binary system, a constant pressure of 8.0 MPa is imposed both on the wellhead of production and injection after the production well is occupied by CO 2 (about 3 years).CO 2 is injected into the reservoir at a constant temperature of 15 ∘ C for 20 years.

Simulation Approach
3.1.Governing Equations.The process of fluid flow in subsurface can be divided into two parts, flow in the wellbores and flow in the reservoir formation.Therefore, we employed the wellbore-reservoir coupled simulator T2well [17][18][19].It is an extension of the general multiphase, multicomponent, nonisothermal simulator TOUGH2 V2 [12].The program assigns the wellbore and reservoir as two subdomains, in which flows are controlled by appropriate laws, respectively.In the reservoir, the flow is described as Darcy's law and in wellbore, it is momentum conservation.To consider a comprehensive description of the thermophysical properties of H 2 O-CO 2 mixtures, the ECO2N V2.0 module was introduced, which reproduces fluid properties including density, viscosity, and specific enthalpy, largely within experimental errors under the temperature and pressure conditions of 10 ∘ C <  < 300 ∘ C,  < 600 bar [20].These fundamental flow equations, used in the T2Well code [18], are summarized in Table 4.In Table 4,  means the accumulation term, and  represents the flux term of mass or energy. represents phase and  is the index for fluid species.For mass conservation,   ,   ,   are density, volumetric fraction, and velocity of phase-, respectively.   is the mass fraction of component  in phase-.For energy conservation, , ℎ, and  are heat conductivity, specific enthalpy of fluid, and inclination angle of wellbore, respectively.For momentum equation, Γ is perimeter of wellbore, and   is the wall shear stress.[21,22] is introduced to calculate the two-phase velocities of CO 2 -water mixtures in wellbores by the following equations.

Drift-Flux Model. Drift-Flux Model
First, the velocity of gas phase can be described by the constitutive relation as below: where the profile parameter  0 is used to account for the effect of local gas saturation and velocity profiles over the pipe (wellbore) cross-section [23];  is the volumetric flux of total mixture, which is defined as Thus, the velocity of liquid could be determined as Then the mixture velocity can be calculated in light of By inserting (3) and ( 5) into ( 6), the volumetric flux  could be described as a function of mixture velocity   and the drift velocity   , as (7) shows: where  #  =    0   + (1 −    0 )  , and it is the profileadjusted average density.The major task is now to calculate the mixture velocity and the drift velocity.

Results and Discussion
4.1.Migration of CO 2 in the Reservoir.The distribution of CO 2 saturation in the reservoir after 3 years of injection is shown in Figure 3.By comparing Figure 3(a) through 3(f) (6 heterogeneous cases) with the homogeneous cases of Figure 3(g), it can be found that the migration of CO 2 is significantly affected by the media (permeability and porosity) heterogeneity.In the homogenous reservoir, the CO 2 saturation distributes circularly around the injection well.However, in the heterogeneous cases, the CO 2 tends to flow through the highly permeable zone and leads to a preferential flow.It bypasses the low permeability zone to reach the production well, even if the flow path is longer (shown in Figures 3(e) and 3(f)).The positions holding high CO 2 saturation match well with highly permeable zone.For the cases that the highly permeable zones locate far away from the production well (Figures 3(c) and 3(d)), there would be more injected CO 2 stay in reservoir instead of being extracted from the production well.

CO 2 Saturation in Production Fluids.
Being injected through the injection wellbore, CO 2 displaced the original water in the reservoir formation, and then water was first produced from the production well.At early stage, the production fluid was pure water, and CO 2 plume expanded within reservoir over time.When it reached the production well, the output fluids became into a mixture of water and CO 2 .As can be seen in Figure 4(b), the time of CO 2 breakthrough was quite different in different cases.In homogeneous case, it takes about 1.3 years for CO 2 to reach the production well, while the breakthrough time of most heterogeneous cases is earlier; for example, it only takes 0.8 years for CO 2 to get through the reservoir and reach the production well in hh2 case, within which there is a highly  permeability zone between two wells.Correspondingly, in hl1 and hl2 cases, the breakthrough time is 1.4 years and 1.5 years, respectively, due to the low permeable belt between the injection and production wells.After CO 2 reached the production well, the saturation of CO 2 in the output fluids increased rapidly to over 0.9 (Figure 4(a)).The dive of curves in the figures is caused by the increasing of production pressure (from 0.1 MPa to 8 MPa as mentioned before).After that, the saturation of CO 2 recovered smoothly and stayed at about 0.92 during the rest operation time.

Injection and Production Rate of Fluids.
The evolution of injection rate of CO 2 of different cases is presented in Figure 5.As shown in the figure, the injection rate of CO 2 rises rapidly after it reaches the production well.When injection and production wellhead pressures are imposed to be 8 MPa, the cyclic pressure difference ( inj −  pro ) reduces to zero from 11.9 MPa, the injection rate of CO 2 instantly fell and stayed relatively constant during the rest time of simulation.The injection rate of CO 2 in hh1 and hh2 is relatively higher compared to other cases due to its high permeability and porosity zone between the two wellbores.Figure 6(b) shows the evolution of production rate of water.It can be seen that the rate increases dramatically like an eruption, when CO 2 appears from the production wellhead.This phenomenon is caused by the CO 2 accumulation in production; when it exceeds a certain quantity, an eruption event would take place.This phenomenon can likely account for some cold geysers eruptions [24,25].After the eruption, the content of water in the mixture fluids goes down.At the instant of raising the production pressure and reducing the injection pressure to 8 MPa, both the production rate of water and CO 2 fall off.Note that the production rate of water declines continually but that of CO 2 rises gradually.Even though the heat transmission fluid is CO 2 , water still takes a large proportion of the production fluids because of the huge density difference in production well.Simulation results show that low permeability zone between the wellbores (hl1 and hl2) does not correspond to the lowest production rate of fluids as we expected.On the contrary, water production rate is relatively larger compared to other cases.This is because the resident water saturation at production well bottom in hl1 and hl2 is higher than that of other cases.On the other hand, the highly permeable belts connecting two wells may cause the high production rate of CO 2 and low production rate of water (hh1 and hh2).As can be seen in the graph, the fluids production rates of hp1 and hp2 are really low during the simulation time, which could affect the economic feasibility of the CPG system.

4.4.
Temperature of Production Fluids.The temperature of production fluids is affected by the flow path and heat exchange within the reservoir.As shown in Figure 7, the temperature of production water distributes between 100 and 115 ∘ C before CO 2 reaches the production wellhead, which is somewhat less than the reservoir initial temperature (120 ∘ C).This temperature drop in the production well is due to the heat exchange (loss) with the surrounding rock around the wellbore.Consequently, when CO 2 begins to produce in the wellhead, the temperature of mixture output fluids drops rapidly.It is caused by the Joule-Thomson expansion in the wellbore with the pressure decrease from the downhole to the wellhead.It can be seen that the temperature drop of output fluids from the downhole to wellhead is approximate 40 ∘ C.  Enhancing the production pressure increases the temperature of output fluids by 10 to 14 ∘ C.
Comparing the temperature of production fluids between different cases, it can be found that, to a certain extent, the temperature of output mixture fluids is proportional to the water output rate.It can be explained from two aspects.Firstly, water can counteract the Joule-Thompson effect and help keep the temperature.Secondly, higher water output rate implicates that more space of reservoir is occupied by CO 2 , and more energy is extracted by the output fluid.That is just the reason why the cases with low permeability zone between the two wells correspond to a higher output fluids temperature (hl1 and hl2).

CO 2 Storage.
To investigate the CO 2 storage capacity in the reservoir, we define the storage rate of CO 2 as the injection rate minus production rate of CO 2 (CO 2inj − CO 2pro ).The storage rate of CO 2 of different cases is shown in Figure 8(a).More CO 2 is stored in the cases with highly permeable zone located away from production well.The storage amount of CO 2 , that is, the storage rate of CO 2 integral of time, is shown in Figure 8(b).Similarly, the storage amount of CO 2 of cases hp1 and hp2 (with highly permeable zone deviating from production well) is relatively higher compared to other cases.However, the storage amount of CO 2 is also affected by cumulative injection of CO 2 .So the storage ratio of CO 2 is calculated by storage amount of CO 2 divided by cumulative injection of CO 2 and is shown in Figure 9.In cases hh1 and hh2, CO 2 tends to migrate to production wellbore through high permeability zone between the two wellbores and the storage ratio of CO 2 in the reservoir is less than homogeneous case.On the contrary, in cases hp1 and hp2 (with high permeability belt deviating from producer) a portion of CO 2 will transport away from production wellbore and trapped in the reservoir.The storage ratio of cases hp1 and hp2 is 52.7% and 55.3%, respectively, which is relatively higher compared to homogeneous case (45%).4.6.Temperature Distribution of Reservoir.The temperature drop in reservoir is caused by heat transfer between rock and fluids.After 20 years of production, the temperature of the rock matrix between the injector and producer is significantly reduced.It is highly affected by the hydrological heterogeneity.In cases hl1 and hl2 (Figures 10(e) and 10(f)), the temperature drop around the production wellbore is less than other cases due to the low permeability and porosity zone between the two wellbores.

The System Heat Extraction Rate.
The net heat extraction rate is calculated by the following [26]: where  is the net heat extraction rate,  is the production flow rate, ℎ inj is the specific enthalpy of injection fluids, and ℎ pro is the specific enthalpy of production fluids.
The simulated heat extraction results (Figure 11) indicate that media heterogeneity affects the heat extraction rate greatly.In cases with low permeability zone between production and injection wells (hl1 and hl2), the heat extraction rates are similar with the homogeneity case and stay relatively constant in the middle and later periods.Cases hp1 and hp2, with high permeability zones extending far away from production well, have low heat extraction rates during the entire simulation time.The other heterogeneous cases with high permeability zone between two wells (hh1 and hh2) obtain the maximum extraction rates at early time (2 to 3 years) and decrease quickly in the later periods.

Concluding Remarks
We have built a three-dimensional wellbore-reservoir coupled model with consideration of permeability and porosity heterogeneity based on the geological and thermal-physical conditions of Songliao Basin, China.A total of 7 case simulations were performed.The following conclusions can be drawn.
Heterogeneity of reservoir's hydrological properties (in terms of permeability and porosity) affects the migration of

Figure 2 :
Figure 2: Permeability distribution of six heterogeneous cases.

Figure 3 :
Figure 3: Carbon dioxide saturation distribution in the cross-section of subdomain in the depth of 3050 m for different cases after 3 years.

Figure 4 :Figure 5 :
Figure 4: Variation of CO 2 saturation in output fluids versus time at the wellhead.

Figure 6 :
Figure 6: Variation of production rate of CO 2 (a) and water (b) versus time for different cases.

Figure 8 :
Figure 8: Storage rate (a) and storage amount (b) of CO 2 within reservoir versus time.

Figure 9 :
Figure 9: Storage ratio of CO 2 in reservoir after 20 years for different cases.

Figure 11 :
Figure 11: Variation of heat extraction rate over time for different cases.

Table 3 :
Parameter setting of different modeling scheme.