PROCEEDINGS, Fortieth Workshop on Geothermal Reservoir Engineering Stanford University, Stanford, California, January 26-28, 2015 SGP-TR-204

Validation of a Numerical Reservoir Model of Sedimentary Geothermal Systems Using Analytical Models JaeKyoung Cho1, Chad Augustine2 and Luis E. Zerpa1 1

Petroleum Engineering Department, Colorado School of Mines, Golden, CO, USA 2 National Renewable Energy Laboratory, Golden, CO, USA [email protected], [email protected], [email protected]

Keywords: sedimentary geothermal, numerical modeling, validation, well doublet, reservoir modeling. ABSTRACT The commercial production of geothermal energy from sedimentary rock formations with relatively low permeability could expand the current geothermal energy resources toward new regions in the U.S. The feasibility of commercial geothermal production from sedimentary reservoirs can be studied using numerical reservoir modeling, which allows studying different well configurations and productivity enhancement techniques. In this paper we present the validation of a numerical reservoir model with respect to an analytical model, and the process followed to achieve an acceptable match between the numerical and analytical solutions. The analytical model used is based on the work of Gringarten (1978), which consists of a conceptual sedimentary geothermal reservoir model, considering an injection and production well doublet in a homogeneous reservoir. The numerical modeling is conducted using a commercial thermal reservoir simulator. In order to reproduce the analytical model results, the numerical simulation model is modified to include the same assumptions of the analytical model. The following model parameters are considered to obtain an acceptable match between the numerical and analytical solutions: grid block size, time step and reservoir areal dimensions; the latter related to boundary effects on the numerical solution. Here we present lessons learnt and guidelines to improve the numerical simulation solutions. The agreement between the numerical and analytical models allows us to proceed with confidence to study reservoir enhancement techniques using numerical reservoir simulation. The validated numerical reservoir model is used to compare the performance of the geothermal system with and without enhancement treatments (e.g., hydraulic fracturing). 1. INTRODUCTION Sedimentary enhanced geothermal systems (SEGS) represent a great potential as an energy resource. Unlike conventional convectiondominated high-enthalpy systems, SEGS are not limited to volcanic or tectonic areas that are less extensive than sedimentary basins in their areal dimension (Bundschuh and Suarez, 2010). In a report led by the Massachusetts Institute of Technology, it was estimated that the US could extract 200,000×1018 Joules out of an EGS total potential of 13,000,000×10 18 Joules for energy utilization (Tester et al., 2006). Sedimentary formations with high temperature at depths ranging from 2 to 6 km can be potential candidates for commercial enhanced geothermal energy production. The permeability of economically viable potential target formations at such depths is expected to be in the range of 1 md to 100 md. There is a wide range of variation in petrophysical properties depending on geologic settings. In low permeability formations ( 10%) to be considered acceptable.

Table 4: Summary of numerical results of thermal breakthrough time obtained with different grid block sizes, and comparison with analytical solution. Grid block size (m)

Reservoir size (m)

Well spacing (m)

100 x 100 50 x 50 10 x 10 5x5 2.5 x 2.5

4,500 x 3,000 4,500 x 3,000 4,500 x 3,000 4,500 x 3,000 4,500 x 3,000

1,500 1,500 1,500 1,500 1,500

Thermal breakthrough time Analytical Numerical solution solution 30.0 16.0 30.0 19.0 30.0 22.0 30.0 23.0 30.0 24.0 5

Relative error (%) 46.7 36.6 26.6 23.3 20.0

Cho, Augustine and Zerpa 3.3. Effect of reservoir no-flow boundaries on thermal breakthrough time In this section, we improve further the numerical solution by considering the effect of the distance of no-flow reservoir boundaries from the wells. Figure 5 shows a comparison between the base reservoir model, with reservoir dimensions used in the previous section, and an enlarged reservoir model, where the no-flow boundaries are further away from the wells. Figure 6 shows the results obtained with a grid block size of 50 x 50 m and different reservoir dimensions. The reservoir dimensions are specified in terms of the distance between the wells. Table 5 summarizes the numerical results obtained with different reservoir model dimensions for a grid block size of 50 x 50 m, and presents a comparison with the analytical solution. The thermal breakthrough time increases with increasing reservoir area, and the relative error between the numerical and analytical solution decreases. By increasing the reservoir model size, the effect of no-flow boundaries on the numerical solution is decreased, and the numerical solution better approximates the analytical. There is a reservoir model area where further increase does not improve the solution significantly.

Figure 5: Base reservoir model with grid block size 50 x 50 m (left: 4500 x 3000 m) and enlarged model (right: 7500 x 6000 m) to study boundary effects on numerical solution.

Temperature (ºC)

160 155 150 145 140 135 0

4500  3000 7550  6050 11550  10050 27550  26050 107550  106050 10

20

Time (yr)

30

40

Figure 6: Plot of temperature at the production well location as function of time, showing effect of reservoir model areal dimension (in meters) on calculated thermal breakthrough time, with a grid block size of 50 x 50 m. Table 5: Summary of numerical results of thermal breakthrough time obtained with different reservoir model dimensions for a grid block size of 50 x 50 m, and comparison with analytical solution.

Reservoir size (m)

Dimension ratio to well spacing (x:y)

4,500 by 3,000 7,500 by 6,000 11,500 by 10,000 27,500 by 26,000 10,7500 by 10,600

3:2 5:4 7.7:6.7 18.3:17.3 71.6:70.6

Thermal breakthrough time (years) Analytical Numerical solution solution 30 19 30 21 30 22 30 22 30 22

Relative error (%) 36.6 30.0 26.6 26.6 26.6

Additionally, we studied the effect of reservoir boundaries on the numerical solution using different grid block sizes. Figure 7 shows the temperature at the production well location as a function of time for selected cases to illustrate the improvement in the numerical solution by extending the reservoir dimension and refining the grid block sizes. Table 6 presents a summary of the numerical results and a comparison with the analytical solution. The match between the numerical and analytical solutions improves with the increase in the 6

Cho, Augustine and Zerpa reservoir model dimensions and the decrease of the grid block size used. An acceptable lowest relative error of 6.67 % is obtained with an extended reservoir model with grid block size of 2.5 x 2.5 m.

Temperature (ºC)

160 155 150 145 140 135 0

10  10 (4500  3000) 10  10 (24510  23010) 2.5  2.5 (37765  36265) 10

20

Time (yr)

30

40

Figure 7: Plot of temperature at the production well location as function of time, showing effect of no-flow reservoir boundaries (reservoir model areal dimension in meters) on calculated thermal breakthrough time, for selected cases with different grid block sizes. Table 6: Summary of numerical results of thermal breakthrough time obtained with different reservoir model dimensions for selected cases with different grid block sizes, and comparison with analytical solution. Grid block size (m)

Reservoir size (m)

10 x 10 10 x 10 2.5 x 2.5

4,500 x 3,000 24,510 x 23,010 37,765 x 36,265

Thermal breakthrough time (years) Analytical Numerical solution solution 30 23 30 26 30 28

Relative error (%) 23.33 13.33 6.67

3.4. Optimum grid configuration The next step was the study of different numerical grid configurations, in order to obtain an optimal grid system with a minimum number of grid blocks and practical computational time. The numerical grid configuration is defined by specifying the grid block size of the interest area in the center of the reservoir model containing the pair of wells, the total reservoir areal dimensions, and the grid block size distribution in the reservoir boundary surrounding the interest area. Table 7 presents the parameters of different numerical grid configurations considered in this study, and Figure 8 shows an areal view of the grid configurations. The optimal numerical grid configuration is selected from the numerical solution that shows the longest thermal breakthrough time. Figure 9 shows the temperature at the production well as a function of time obtained using the different grid configurations. Grid (f) shows the best solution and it is associated with a practical computational time (from Table 7). Table 7. Parameters of numerical grid configurations studied. Grid ID

Grid block size of interest area (m)

Reservoir size (m)

(a) (b)

20 x 20 20 x 20

4,500 x 3,000 34,500 x 33,000

(c)

20 x 20

34,500 x 33,000

(d) (e)

20 x 20 20 x 20

34,500 x 33,000 66376 x 64876

(f)

20 x 20

34,500 x 33,000

(g)

20 x 20

34,500 x 33,000

Grid block size distribution of reservoir boundary N/A Coarser Linear, 20 x 20, for coordinates of interest area Logarithmic coarser Logarithmic finer Logarithmic medium 20 x 20 entire reservoir model

7

Number of grid blocks

Computational time

33,751 33,759

6m 44s 6m 45s

34,504

6m 57s

34,191 42,650

7m 3s 8m 16s

34,591

7m 00s

2,881,725

11h 22m 02s

Cho, Augustine and Zerpa Applying the optimal numerical grid configuration, grid (f), to several cases with different grid block sizes in the interest area, we verify that the calculated thermal behavior of the geothermal reservoir is more accurate when using a finer grid. Figure 10 shows the temperature change at the production well with time for different grid block sizes using grid (f). Similar to Figure 4, the solution with the optimal grid configuration also shows longer thermal breakthrough times with smaller grid block sizes, and a sharper thermal front arriving to the production well, which indicates a decrease of the numerical dispersion. Table 8 presents a summary of the results obtained with different grid block sizes and the optimal grid configuration. With the two smallest grid block sizes (5 x 5 m, and 2.5 x 2.5 m) the relative error with respect to the analytical solution is reduced to 6.6%.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 8. Grid configurations considered in this study.

Temperature (ºC)

160 155 150 145 140 135 0

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

10

20

Time (yr)

30

40

Figure 9. Plot of temperature at the production well location as function of time, for different numerical grid configurations.

8

Cho, Augustine and Zerpa

Temperature (ºC)

160 155 150 145

100  100 50  50 20  20 10  10 55 22

140 135 0

10

20

Time (yr)

30

40

Figure 10. Plot of temperature at the production well location as function of time, for different grid block sizes using grid configuration (f).

Table 8. Summary of numerical results of thermal breakthrough time obtained with different grid block sizes using grid (f), and comparison with analytical solution. Grid block size (m)

Reservoir size (m)

Well spacing (m)

100 x 100 50 x 50 20 x 20 10 x 10 5x5 2x2

34,500 x 33,000 34,500 x 33,000 34,500 x 33,000 34,500 x 33,000 34,500 x 33,000 34,500 x 33,000

1,500 1,500 1,500 1,500 1,500 1,500

Thermal breakthrough time Analytical Numerical solution solution 30 19 30 22 30 25 30 27 30 28 30 28

Relative error (%) 36.6% 26.6% 16.6% 10.0% 6.6% 6.6%

3.5. Time step effect We perform an evaluation of the time step used during the numerical solution of the flow and thermal behavior of the geothermal reservoir, which also plays a role in the accuracy of the solution by affecting the numerical dispersion. Figure 11 shows the temperature change at the production well with time obtained using different time step values for a grid block size of 10 x 10 m. For time steps less than 1 month the change in the solution is not significant, which indicates that a time step of 1 month is an appropriate value for the maximum time step to be used in the simulation.

Temperature (ºC)

160 155 150 T T T T

145 140 135 0

10

= 10 days = 30 days = 365 days = 1000 days 20

Time (yr)

30

40

Figure 11. Plot of temperature at the production well location as function of time, for different time steps used during the numerical simulation with a grid block size of 10 x 10 m.

9

Cho, Augustine and Zerpa 4. STUDY OF RESERVOIR ENHANCEMENT TECHNIQUES For the evaluation of reservoir enhancement techniques using the numerical reservoir model, the initial assumptions made for the comparison with analytical model are modified to account for a more realistic scenario. The following assumptions are used:        

Horizontal aquifer with uniform thickness Sealed aquifer by cap/bed rocks Constant injection rate equal to production rate (i.e., steady-state condition) Initially, rock and fluid are in thermal equilibrium at same temperature Horizontal heat conduction is considered Vertical heat conduction from cap/bed rock is considered Water properties are function of pressure and temperature Thermal conductivity of water and rock are constant

The base case consists of a well doublet system with a vertical injection well and vertical production well with open hole completions through the entire thickness of the target sedimentary reservoir and separated by a distance (well spacing) of 1,500 m. 4.1. Reservoir enhancement techniques Sedimentary geothermal systems with reservoir enhancement techniques were modeled and their performance compare to that of the base case described above. The following well configurations were implemented to represent the reservoir enhancement techniques considered in this study: 1. 2. 3. 4.

Vertical wells doublet with hydraulic fractures (shown in Figure 12a). Hydraulic fractures half-length of 500 m, and fractures height of 50 m. Horizontal wells with open-hole completion (shown in Figure 12b). Length of horizontal wells is 1000 m. Horizontal well is located in center of formation thickness. Horizontal wells with longitudinal hydraulic fracture (shown in Figure 12c). Length of horizontal wells is 1000 m, and height of longitudinal fracture is 50 m along the length of the horizontal well (25 m fracture half-length). Horizontal wells with multi-stage hydraulic fractures (shown in Figure 12d). Length of horizontal wells is 1000 m, with six evenly spaced hydraulic fractured stages, fractures half-length of 100 m, and fractures height of 50 m (i.e., fractures span entire height of formation). Horizontal well is located in center of formation thickness.

Figure 12: Four different reservoir enhancement techniques considered in this study: a) vertical well doublet with hydraulic fractures; b) horizontal wells with open-hole completions; c) horizontal wells with longitudinal hydraulic fractures; and d) horizontal wells with multi-stage hydraulic fractures. 4.2. Comparison of base case and reservoir enhancement techniques The performance of the different scenarios considered is evaluated in terms of the hydraulic behavior (i.e., well productivity/injectivity), and thermal evolution of the reservoir (i.e., thermal breakthrough time). 10

Cho, Augustine and Zerpa 4.2.1. Reservoir hydraulic behavior This section considers the hydraulic behavior of the reservoir in terms of the pressure difference between the injection and production wells, and the well productivity index (a measure of the production or injection capacity of a well). The reservoir hydraulic behavior is improved by the use of hydraulic fractures, and further improved by the use of horizontal wells. Table 9 presents the estimated reservoir hydraulic parameters for the base case and the cases with reservoir enhancement techniques. Since all cases consider the injection rate equal to the production rate, with a reservoir voidage factor approximating unity, the average reservoir pressure approximates the same value in all cases. However, the local reservoir pressure around the injection and production wells does vary among the cases, with the pressure highest around the injection well and lowest around the production well as one would expect. The difference between the injection and production well pressures is largest for the base case. The injection well bottomhole pressure decreases, while the production well bottom-hole pressure increases, with the introduction of reservoir enhancement techniques. This results in the reduction of the pressure drawdown between the wells beyond the base case for the cases with reservoir enhancement techniques. The productivity index for injector and producer wells increases with the use of enhancement techniques, which represents the improvement of the injection/production capacity of the well by increase in the surface area connecting the well and the reservoir rock matrix through the fractures with each incremental case. Table 9. Estimated reservoir hydraulic behavior parameters for base case and reservoir enhancement techniques cases.

Well Configuration

Bottomhole Average pressure reservoir injection well pressure (kPa) (kPa)

Bottomhole pressure production well (kPa)

Pressure drawdown, DP (kPa)

Productivity Index Injector (l/s-bar)

Producer (l/s-bar)

Base case – well doublet Vertical + hydraulic fracture

28,374

38,349

22,865

15,484

0.94

1.70

28,364

30,523

26,659

3,864

4.33

5.48

Horizontal open-hole

28,367

30,231

27,095

3,136

5.01

7.34

28,368

29,956

27,241

2,715

5.88

8.29

28,369

29,795

27,325

2,470

6.55

8.94

Horizontal with longitudinal fracture Horizontal with multistage fractures

4.2.2. Reservoir thermal behavior This section considers the thermal behavior of the reservoir in terms of the temperature at the location of the production well, and thermal breakthrough time (time when temperature begins to decrease at production well). Figure 13 shows the temperature change with time at the location of the production well obtained with the base case (vertical well doublet) and the different reservoir enhancement techniques considered. The typical behavior of temperature as a function of time observed in each solution presents a constant temperature equal to the initial reservoir temperature, followed by a decline in temperature because of arrival of colder water to the production well. After this time the temperature at the location of the production well decreases constantly. The thermal breakthrough time increases with the introduction of reservoir enhancement techniques (Table 10). The use of vertical wells with hydraulic fractures improves the thermal behavior of the geothermal reservoir by increasing the thermal breakthrough time relative to the base case. Greater thermal breakthrough times are obtained with the horizontal well cases than the vertical well cases. The use of horizontal wells reduces the velocity of the cold front moving within the reservoir, hence increasing the thermal breakthrough time when compared to vertical well doublet cases. 162

Temperature (ºC)

160 158 156 154 152 150 148 146 0

Vert. doublet Vert. frac. Hor. open-hole Hor. longitudinal frac. Hor. transverse frac. 10

20

30

Time (yr)

40

50

Figure 13. Plot of temperature at the production well location as function of time, for base case and reservoir enhancement techniques cases. 11

Cho, Augustine and Zerpa Table 10. Thermal breakthrough time obtained from base case and different reservoir enhancement techniques. Well Configuration Base case – well doublet Vertical + hydraulic fracture Horizontal open-hole Horizontal with longitudinal fracture Horizontal with multi-fracture

Thermal breakthrough time (yr) 27 36 40 42 41

5. CONCLUSIONS A rigorous study comparing the results of the analytic model for the hydraulic and thermal behavior of flow between a well doublet pair in a sedimentary reservoir to numerical modeling of the same system showed that the numeric model results for the hydraulic behavior are insensitive to the grid configuration, but that the model results for thermal behavior can be heavily dependent on grid size, modeled reservoir area (boundary effects), time step, and grid configuration. Systematic model runs showed that insufficient grid sizing causes numerical dispersion that causes the numerical model to underestimate the thermal breakthrough time compared to the analytic model. As grid sizing is decreased, the model results converge on a solution (Figure 4). Likewise, insufficient reservoir model area introduces boundary effects in the numerical solution that cause the model results to differ from the numerical model. Increasing the reservoir model area improves the agreement of the numerical model results to the analytical model results by better approximating the “infinite reservoir” assumption of the analytical model. As reservoir model area increases, the numerical model results converge and the benefits of using even larger reservoir model areas diminishes (Figure 6). A numerical model using a small grid size (2.5x2.5 m) and large reservoir model area (37,765m x 36,265m) was found to give numerical results for the thermal breakthrough time in good agreement with the analytical model (see Table 6 and Figure 7). However, using small grid sizes and large reservoir model areas results in a large number of grid cells and prohibitive computational times (see Table 7, case g). Different grid configurations that use a fine grid mesh in the area of interest and a coarse grid mesh in the outer regions where the temperature and pressure gradients are small were tested to find one that gave an accurate solution in a reasonable amount of computational time. The resulting optimal grid configuration was found to use a total reservoir model area of 34.5 km x 33 km, a 5x5 m spacing in the area of interest closest to the well doublet, and a logarithmic medium algorithm to determine grid block size distribution in areas around the reservoir boundary. This configuration resulted in a computed temperature breakthrough time within 6.6% of the analytic solution. A numerical simulation using 2x2 m spacing in the area of interest confirmed that additional reductions in grid size did not significantly change the numerical model results. Finally, a study of the impact of time step found that assuming a one-month time step gave consistent results. The study demonstrated that special attention must be paid to the reservoir grid used to ensure that numerical results are accurate and not overly influenced by numerical dispersion and boundary effects. Trade-offs must be made in terms of the number of grid cells used in the model and the computational time required. The systematic study of the influence of these factors on numerical model results compared to an analytic model demonstrated that our numerical model can give accurate results and gave us guidance on the reservoir model area and grid sizing to use in future numerical modeling. Four different reservoir enhancement techniques were compared against a base case. The base case consists of a sedimentary geothermal reservoir with an injection and production well doublet. The performance of the different scenarios considered was evaluated in terms of the hydraulic behavior (i.e., productivity/injectivity), and thermal evolution of the reservoir (i.e., thermal breakthrough time). The reservoir enhancement techniques can have a significant impact on geothermal reservoir performance. The productivity index of injector and producer wells increase by a factor of 5, approximately, with the use of enhancement techniques. This is due to the augmented surface area connecting the well and the reservoir rock matrix through the fractures in each reservoir enhancement technique tested. The geothermal reservoir lifetime is increased by 50%, approximately, by the use of horizontal wells using the same well spacing as the base case (vertical well doublet) and same injection/production rate. The use of horizontal wells reduces the velocity of the cold front moving within the reservoir, hence increasing the thermal breakthrough time when compared to vertical well doublet cases. The most promising reservoir enhancement techniques are horizontal wells with longitudinal fracture and horizontal wells with multi-stage hydraulic fractures. ACKNOWLEDGEMENTS This work was supported by the U.S. Department of Energy, Office of Energy Efficiency and Renewable Energy (EERE), Geothermal Technologies Office (GTO) under Contract No.DE-AC36-08-GO28308 with the National Renewable Energy Laboratory. REFERENCES Allis, R. et al., 2011. The potential for basin-centered geothermal resources in the Great Basin. Geothermal Resources Council Transactions, 35: 683-688. Augustine, C., 2014. Analysis of sedimentary geothermal systems using an analytical reservoir model. Geothermal Resources Council Transactions, 38: 641-647. Bjelm, L. and Alm, P., 2010. Reservoir Cooling After 25 Years of Heat Production in the Lund Geothermal Heat Pump Project, Proceedings, World Geothermal Congress. 12

Cho, Augustine and Zerpa Bundschuh, J. and Suarez, M.C., 2010. Introduction to the numerical modeling of groundwater and geothermal systems: Fundamentals of mass, energy and solute transport in poroelastic rocks. Multiphysics Modeling. CRC Press. Deo, M., Roehner, R., Allis, R. and Moore, J., 2013. Reservoir modeling of geothermal energy production from stratigraphic reservoirs in the Great Basin, Thirty-Eighth Workshop on Geothermal Reservoir Engineering. Stanford University, Stanford, CA. Gringarten, A., 1978. Reservoir lifetime and heat recovery factor in geothermal aquifers used for urban heating. pure and applied geophysics, 117(1-2): 297-308. Gringarten, A.C. and Sauty, J.P., 1975. A theoretical study of heat extraction from aquifers with uniform regional flow. Journal of Geophysical Research, 80(35): 4956-4962. Muskat, M. and Wyckoff, R.D., 1946. The flow of homogeneous fluids through porous media. International series in physics. J. W. Edwards, In., Ann Arbor, Mich., xix, 763 p. pp. Satman, A., 2011. Sustainability of Geothermal Doublets, Proceedings, Thirty-Sixth Workshop on Geothermal Reservoir Engineering, Stanford University, Stanford, California, Jan. Tester, J.W. et al., 2006. The future of geothermal energy - Impact of Enhanced Geothermal systems (EGS) on the US in the 21st century, Idaho National Laboratory, Idaho Falls, Idaho.

13