Investigative Coupled Thermo-Hydro-Mechanical Modelling Approach for Geothermal Heat Extraction through Multistage Hydraulic Fracturing from Hot Geothermal Sedimentary Systems

The meaningful utilization of artificially created multiple fractures in tight formations is associated with the performance behavior of such flow channels, especially in the case of thermal energy extraction from sedimentary geothermal system. In this study, an innovative idea is presented to develop a numerical model for geothermal energy production based on concrete physical performance of an artificially created tensile multi-fracture system in a simplified manner. The state-of-the-art software FLAC3Dplus-TOUGH2MP-TMVOC are integrated to develop a coupled thermo-hydro-mechanical (THM) fictive model for constructing a multi-fracture scheme and estimating heat extraction performance. By incorporating the actual fracture width of newly created subsequent fracture under the effect of stress shadow, cubic law is implemented for fluid flow and geothermal energy production. The results depict that fracture spacing plays a vital role in the energy contribution through multiple fractures. Afterwards, a field case study to design huge multiple hydraulic fractures was performed in the geothermal well GB X1 in North Germany. The attenuation of fracture propagation becomes more significant when massive multiple fracturing operation is performed especially in the case of lower fracture spacing. The fictive model results will be extended to study the geothermal utilization of the North German basin through massive multiple fractures in our future work.


Introduction
The interminable increase in energy demand has triggered the exploitation of unconventional energy resources in recent times. There is an immense need of utilizing renewable energy resources along with dependency on natural fossil fuels. In this recent era of global warming, researchers are more focused to work on efficient, environmentally friendly, and renewable energy resources. Based on these requirements, geothermal energy is among the leading positions, as it provides heat energy which is independent of weather and at the same time is sustainable, renewable, and virtually without greenhouse gas emissions [1][2][3][4][5]. Geothermal energy evolves due to two sources, (a) transfer of energy from hot molten core to exterior of the earth, (b) decay of the radio-active elements [6]. The interior of earth has a temperature above 4000 K (approx.) [7] and is considered to be a huge source of geothermal Although some researchers [33][34][35] have studied and analyzed the geometry and propagation of fractures employing a simultaneous multi-stage hydraulic fracturing technique for economical production of fossil fuels in tight formations, the impact of stress shadow on multiple fracture configuration and consequently for geothermal energy production especially incorporating THM coupling effects still needs critical investigation.

Simulation Concept
The numerical simulation of the hydraulic fracturing process carries out many complex physical processes which occur during mechanical deformation of rock. Mathematical formulation during these processes includes a combination of equation of motion (Equation (1)), continuum equation (Equation (2)), and constitutive equation (Equation (3)) [36].
where q i is heat flow (W/m 2 ) in i direction (x, y, z), λ is thermal conductivity (W/m. • C), T is temperature ( • C), q v is heat source of volume (W/m 3 ), H is stored heat per unit volume (J/m 3 ), c v is specific heat capacity (J/kg. • C), q n component of flux normal to the boundary (W/m 2 ), h is convective heat transfer coefficient (W/m 2 • C), T is the temperature of solid boundary surface ( • C), T e is temperature of surrounding fluid ( • C). In this paper, TOUGH2MP-TMVOC is introduced to simulate the multi-phase, multi-component fluid flow in both fracture and reservoir with the following fluid flow equation (Equation (8)) and mass conservation equation (Equation (9)), respectively.
The permeability of each fractured zone is calculated using cubic law (Equation (10)), as the fluid flow through fracture plane is treated as flow between parallel plates.
where k denotes the fracture permeability (m 2 ), w corresponds to fracture width (m), and f is categorized as fracture roughness value (-). In our model, hydro-mechanical (HM) and thermo-hydro (TH) coupled equations are solved based on numerical computational methods by following the sequential flow scheme as depicted in Figure 1. It is hard to find the work done in the literature, which incorporates the stress shadowing impact, while designing multiple fractures in enhanced geothermal systems. The uniqueness of this study is exemplified by a couple of prominent features. The impact of the previous fracture on each newly created fracture was analyzed based on configured fracture geometry, during multiple hydraulic fracturing jobs, encompassing the stress re-orientation (HM coupling). Afterwards, flow calculations for geothermal utilization were then performed based on cubic law, while integrating the changes that occurred during the fracturing job (TH coupling). While, mechanical changes during the production stage were assumed as negligible. This concept can provide us a good understanding of the preferable flowing paths and a fair estimation of geothermal energy production.
Energies 2020, 13, x FOR PEER REVIEW 4 of 24 In our model, hydro-mechanical (HM) and thermo-hydro (TH) coupled equations are solved based on numerical computational methods by following the sequential flow scheme as depicted in Figure 1. It is hard to find the work done in the literature, which incorporates the stress shadowing impact, while designing multiple fractures in enhanced geothermal systems. The uniqueness of this study is exemplified by a couple of prominent features. The impact of the previous fracture on each newly created fracture was analyzed based on configured fracture geometry, during multiple hydraulic fracturing jobs, encompassing the stress re-orientation (HM coupling). Afterwards, flow calculations for geothermal utilization were then performed based on cubic law, while integrating the changes that occurred during the fracturing job (TH coupling). While, mechanical changes during the production stage were assumed as negligible. This concept can provide us a good understanding of the preferable flowing paths and a fair estimation of geothermal energy production.

Numerical Simulation
In this study, numerical simulator FLAC3D plus was used for hydraulic fracturing based on the method developed by Zhou et al. [38], in which fracture propagates under three-dimensional stress state with fully hydro-mechanical coupling effects among fracture and matrix mediums.

Numerical Simulation
In this study, numerical simulator FLAC3D plus was used for hydraulic fracturing based on the method developed by Zhou et al. [38], in which fracture propagates under three-dimensional stress state with fully hydro-mechanical coupling effects among fracture and matrix mediums.

Fictive Model Generation
To apply the proposed idea, a fictive 1 2 3D model is generated as shown in Figure 2a. A pay zone section (sandstone) of 100 m height is sandwiched between cap and basement rocks having 50 m each, in vertical height/thickness. The model has dimensions of 250 (x) × 450 (y) × 200 m (z) and distributed into 105,000 rectangular elements. Mechanical and hydraulic properties of layers are given in Table 1, while the initial stress and pore pressure conditions are shown in Figure 2b. For the application of multistage fracturing, a horizontal well is located at the center of the pay zone layer at a depth of −3100 m. In addition, the direction of the horizontal well is along the minimum horizontal stress (along the y-axis) to facilitate the development of tensile fractures. With the technique of sequential hydraulic fracturing, four fractures are created into the representative model starting from toe to heel (Figure 3a). A fluid volume of about 650 m 3 is injected into each fracture at an injection rate of 7.2 m 3 /min. The distance between the consecutive fractures is taken as 80 m as a base case.

Fictive Model Generation
To apply the proposed idea, a fictive ½ 3D model is generated as shown in Figure 2 a. A pay zone section (sandstone) of 100 m height is sandwiched between cap and basement rocks having 50 m each, in vertical height/thickness. The model has dimensions of 250 (x) × 450 (y) × 200 m (z) and distributed into 105,000 rectangular elements. Mechanical and hydraulic properties of layers are given in Table 1, while the initial stress and pore pressure conditions are shown in Figure 2 b. For the application of multistage fracturing, a horizontal well is located at the center of the pay zone layer at a depth of −3100 m. In addition, the direction of the horizontal well is along the minimum horizontal stress (along the y-axis) to facilitate the development of tensile fractures. With the technique of sequential hydraulic fracturing, four fractures are created into the representative model starting from toe to heel (Figure 3 a). A fluid volume of about 650 m 3 is injected into each fracture at an injection rate of 7.2 m 3 /min. The distance between the consecutive fractures is taken as 80 m as a base case.     Figure 3b shows the four created fractures widths while minimum horizontal stress shadow is presented in Figure 3c. It is interesting to observe that each fracture has a different fracture width despite the same injection strategy, while half lengths of each generated fracture are the same. In addition, the horizontal stress decreases perpendicular to the fracture plane and interacts with the stress shadow of the next succeeding fracture. Although, stress re-orientation affected the fracture widths of subsequent fractures, but it is observed that all four fractures have the same heights, which also proves a good agreement of strong upper and lower barrier conditions.
In order to further investigate the stress re-orientation effects on fracture width, two simulations are conducted with 60 and 120 m fracture spacing using the same injection strategy as of the base case (80 m). Comparative results of fracture width and minimum horizontal stress shadow are presented in Figure 4. It is noticed that the fracture width changes radically in the case of 60 m spacing as compared to 120 m spacing, while having the same half fracture lengths. Hence, stress shadow has become a source of fracture width variation in a multiple fracture system.   Figure 3 c. It is interesting to observe that each fracture has a different fracture width despite the same injection strategy, while half lengths of each generated fracture are the same. In addition, the horizontal stress decreases perpendicular to the fracture plane and interacts with the stress shadow of the next succeeding fracture. Although, stress re-orientation affected the fracture widths of subsequent fractures, but it is observed that all four fractures have the same heights, which also proves a good agreement of strong upper and lower barrier conditions.
In order to further investigate the stress re-orientation effects on fracture width, two simulations are conducted with 60 and 120 m fracture spacing using the same injection strategy as of the base case (80 m). Comparative results of fracture width and minimum horizontal stress shadow are presented in Figure 4. It is noticed that the fracture width changes radically in the case of 60 m spacing as compared to 120 m spacing, while having the same half fracture lengths. Hence, stress shadow has become a source of fracture width variation in a multiple fracture system.
According to cubic law, the permeability of the created fracture is directly dependent on fracture width. With the variation in fracture width, the flow conductivity of the fracture cannot be the same. However, through fictive modeling results, the relation between stress shadow and fracture conductivity has been explained well, which will be presented in the next section.

Equivalent Model for Geothermal Exploitation
After fracturing operations, production modeling is analyzed further. While, considering the generated fractures as porous and permeable medium, four fracture zones are generated using an equivalent fracture concept based on fracturing results with FLAC3D Plus and then the model is imported to TOUGH2MP-TMVOC for hydro-thermal simulation (as explained in Figure 1). This approach provides benefits of accounting for coupled thermo-hydro-mechanical effects in a time efficient manner. Moreover, TOUGH2MP-TMVOC provides the ability of multiphase multicomponent fluid flow through fractures and reservoir formations, facilitating the temperature and pressure dependent enthalpy, as well [39]. According to cubic law, the permeability of the created fracture is directly dependent on fracture width. With the variation in fracture width, the flow conductivity of the fracture cannot be the same. However, through fictive modeling results, the relation between stress shadow and fracture conductivity has been explained well, which will be presented in the next section.

Equivalent Model for Geothermal Exploitation
After fracturing operations, production modeling is analyzed further. While, considering the generated fractures as porous and permeable medium, four fracture zones are generated using an equivalent fracture concept based on fracturing results with FLAC3D Plus and then the model is imported to TOUGH2MP-TMVOC for hydro-thermal simulation (as explained in Figure 1). This approach provides benefits of accounting for coupled thermo-hydro-mechanical effects in a time efficient manner. Moreover, TOUGH2MP-TMVOC provides the ability of multiphase multi-component fluid flow through fractures and reservoir formations, facilitating the temperature and pressure dependent enthalpy, as well [39].
In this study, a single production well is passed through the created fractures (see Figure 5), instead of using two production wells. The single production well drilled at the center of fracture can provide the energy from both sides of the created fracture. This doublet heat exchanger concept through a single production well will eliminate the expenditures of the second production well. The production points are at a depth of 3060 m, which is about 40 m above the injection points. The temperature from top to bottom of the model varies from 204 to 213 °C, having a geothermal gradient of 0.05 °C/m with surface temperature adjustment. Furthermore, in order to have good fracturing conditions in the flow model, nearby zones are created with smooth permeability decrease in x, y, and z directions. The sensitivity analysis of the nearby zone's permeability change on fluid flow behavior will be discussed in our future work.

Spatial Distribution of Temperature Decline Results During Production
An injection rate of 150 l/min (injected fluid temperature is 55 °C) is used at every injection point for a period of 20 years and correspondingly production rate and enthalpy, which fundamentally depends on the flow conductivity of individual fracture, is measured. Figure 6 shows the temporal evolution of temperature decline in the 1st fracture plane in the x-z direction after a period of 1, 5, 15, and 20 years of production. Similar temperature decline trends are also observed for other fracture planes as well. While, comparing the results during 20 years of production, initially the lower portion of the fracture depletes, and the upper part of the formation depletes afterward. This is due to the fluid movement under the gravitational effect and higher cold water density [1]. The arrangement of injection at lower points, while production from upper points provides an added advantage for harnessing energy from the lower hot section of the model. The temperature from top to bottom of the model varies from 204 to 213 • C, having a geothermal gradient of 0.05 • C/m with surface temperature adjustment. Furthermore, in order to have good fracturing conditions in the flow model, nearby zones are created with smooth permeability decrease in x, y, and z directions. The sensitivity analysis of the nearby zone's permeability change on fluid flow behavior will be discussed in our future work.

Spatial Distribution of Temperature Decline Results during Production
An injection rate of 150 L/min (injected fluid temperature is 55 • C) is used at every injection point for a period of 20 years and correspondingly production rate and enthalpy, which fundamentally depends on the flow conductivity of individual fracture, is measured. Figure 6 shows the temporal evolution of temperature decline in the 1st fracture plane in the x-z direction after a period of 1, 5, 15, and 20 years of production. Similar temperature decline trends are also observed for other fracture planes as well. While, comparing the results during 20 years of production, initially the lower portion of the fracture depletes, and the upper part of the formation depletes afterward. This is due to the fluid movement under the gravitational effect and higher cold water density [1]. The arrangement of injection at lower points, while production from upper points provides an added advantage for harnessing energy from the lower hot section of the model. A comparison of the spatial evolution of temperature decline in the developed model having different fracture spacing is shown in Figure 7. It is observed that, with the passage of time, more area especially in the y-direction depletes significantly and contributes to heat conduction during fluid flow through fractures. Furthermore, with the increase in fracture spacing, a larger area is influenced by cold water injection and hence more stimulated volume in terms of heat conduction is achieved.  A comparison of the spatial evolution of temperature decline in the developed model having different fracture spacing is shown in Figure 7. It is observed that, with the passage of time, more area especially in the y-direction depletes significantly and contributes to heat conduction during fluid flow through fractures. Furthermore, with the increase in fracture spacing, a larger area is influenced by cold water injection and hence more stimulated volume in terms of heat conduction is achieved. A comparison of the spatial evolution of temperature decline in the developed model having different fracture spacing is shown in Figure 7. It is observed that, with the passage of time, more area especially in the y-direction depletes significantly and contributes to heat conduction during fluid flow through fractures. Furthermore, with the increase in fracture spacing, a larger area is influenced by cold water injection and hence more stimulated volume in terms of heat conduction is achieved.

Heat Production Results and Discussion
The difference in injection and production enthalpies provides a good measure of heat production through the simulated model. The produced net-energy can be calculated using the following relationship [40]: where q is the production rate (kg/s), h i is the enthalpy of injected water (j/kg), and h o is the produced fluid enthalpy (j/kg). Figure 8 shows the comparison of temperature decline trends within 20 years of the time period in four individual fractures with different fracture spacing. The short distance between the injection and production point causes an early temperature breakthrough and a sharp decrease in temperature trend. It is quite evident that the breakthrough at the production end can be delayed with the increase in distance of injection and production location.

Heat Production Results and Discussion
The difference in injection and production enthalpies provides a good measure of heat production through the simulated model. The produced net-energy can be calculated using the following relationship [40]: where q is the production rate (kg/s), hi is the enthalpy of injected water (j/kg), and ho is the produced fluid enthalpy (j/kg). Figure 8 shows the comparison of temperature decline trends within 20 years of the time period in four individual fractures with different fracture spacing. The short distance between the injection and production point causes an early temperature breakthrough and a sharp decrease in temperature trend. It is quite evident that the breakthrough at the production end can be delayed with the increase in distance of injection and production location. The corresponding comparison of net-energy contribution (%) for each fracture after 1, 5, and 20 years of production with variable spacing is illustrated in Figure 9. It can be seen that, with 60 m spacing, fracture 1 is contributing least, while fracture 4 is contributing maximum in terms of energy production. This trend remains the same throughout the 20 years of life span. In the case with 80 m fracture spacing, the contribution rate somewhat becomes homogenous as compared to the 60 m spacing case. However, fractures 1 and 2 are still contributing less compared to fractures 3 and 4. An analogous heat contribution is observed with 120 m spacing, which identifies enough spacing for homogeneous heat production through each fracture. It is observed that the unevenness of energy contribution through each fracture is amplified due to shorter distance among fractures in the case of 60 m spacing. The corresponding comparison of net-energy contribution (%) for each fracture after 1, 5, and 20 years of production with variable spacing is illustrated in Figure 9. It can be seen that, with 60 m spacing, fracture 1 is contributing least, while fracture 4 is contributing maximum in terms of energy production. This trend remains the same throughout the 20 years of life span. In the case with 80 m fracture spacing, the contribution rate somewhat becomes homogenous as compared to the 60 m spacing case. However, fractures 1 and 2 are still contributing less compared to fractures 3 and 4. An analogous heat contribution is observed with 120 m spacing, which identifies enough spacing for homogeneous heat production through each fracture. It is observed that the unevenness of energy contribution through each fracture is amplified due to shorter distance among fractures in the case of 60 m spacing.
The rise in temperature of fluid depends on the time duration, which it spends in fracture before production. If the fluid stays in fracture for a longer period, the temperature of the fluid rises more due to heat conduction through the hot surrounding environment. Higher fracture permeability provides a chance for fluid to flow further away from the injection point parallel to the fracture plane before reaching to the producing end, absorbing more reservoir heat. The comparative results of cumulative energy produced and heat power production with time through individual fractures with different spacing is presented in Figure 10. In the case of 60 m spacing, fluid is produced too early without spending much time in fractures 1 and 2 as compared to fractures 3 and 4. The rise in temperature of fluid depends on the time duration, which it spends in fracture before production. If the fluid stays in fracture for a longer period, the temperature of the fluid rises more due to heat conduction through the hot surrounding environment. Higher fracture permeability provides a chance for fluid to flow further away from the injection point parallel to the fracture plane before reaching to the producing end, absorbing more reservoir heat. The comparative results of cumulative energy produced and heat power production with time through individual fractures with different spacing is presented in Figure 10. In the case of 60 m spacing, fluid is produced too early without spending much time in fractures 1 and 2 as compared to fractures 3 and 4. Energy tends to become closer for each fracture in the case of 80 m distance, while 120 m spacing is quite enough to compensate for the fracture permeability variations in terms of less interference among the fracture fluid flow area to contribute to heat production. Figure 11 and Figure 12 shows the cumulative energy produced and total production power from combined four half fracture areas during 20 years of production with different fracture spacing, respectively. From the beginning, maximum energy is produced with 120 m fracture spacing, which corresponds to maximum thermal production power of 3.7 MW as compared to 80 and 60 m spacing, while the difference in this contribution rate becomes less till the end of 20 years of production. Energy tends to become closer for each fracture in the case of 80 m distance, while 120 m spacing is quite enough to compensate for the fracture permeability variations in terms of less interference among the fracture fluid flow area to contribute to heat production. Figures 11 and 12 shows the cumulative energy produced and total production power from combined four half fracture areas during 20 years of production with different fracture spacing, respectively. From the beginning, maximum energy is produced with 120 m fracture spacing, which corresponds to maximum thermal production power of 3.7 MW as compared to 80 and 60 m spacing, while the difference in this contribution rate becomes less till the end of 20 years of production. Energy tends to become closer for each fracture in the case of 80 m distance, while 120 m spacing is quite enough to compensate for the fracture permeability variations in terms of less interference among the fracture fluid flow area to contribute to heat production. Figure 11 and Figure 12 shows the cumulative energy produced and total production power from combined four half fracture areas during 20 years of production with different fracture spacing, respectively. From the beginning, maximum energy is produced with 120 m fracture spacing, which corresponds to maximum thermal production power of 3.7 MW as compared to 80 and 60 m spacing, while the difference in this contribution rate becomes less till the end of 20 years of production. The fictive model results having small fracture areas may not be similar in the case of huge fracturing operations but surely have provided good grounds to analyze large geothermal reservoirs considering stress shadow impact. In continuation, a case study was performed with huge multiple fracturing operations using the field data of the North German basin in the subsequent section.

Field Case Study
The inception of the GB X1 project Germany deals with the progressive results obtained from the test well in the nearby area which has provided basics for the second well named GB X1. The purpose of this project was the utilization of high temperature tight sedimentary formations for geothermal energy production of about 2 MW thermal to fulfill the energy requirement of surrounding residents [41]. Massive single hydraulic fracturing operation was performed on the targeted formation of Buntsandstein at a depth of 3660 m, by injecting 20,000 m 3 of fresh water carrying no proppants to acquire the fracture area of 1.1 km 2 . The whole fracturing operation was performed in around 5 days with several injection pauses. Several tests were also performed to check stress levels and fracture development and found that the fracture had retained its high conductivity without any proppant usage [41,42].
In our case study, a single hydraulic fracture was used for the verification of the simulated model The fictive model results having small fracture areas may not be similar in the case of huge fracturing operations but surely have provided good grounds to analyze large geothermal reservoirs considering stress shadow impact. In continuation, a case study was performed with huge multiple fracturing operations using the field data of the North German basin in the subsequent section.

Field Case Study
The inception of the GB X1 project Germany deals with the progressive results obtained from the test well in the nearby area which has provided basics for the second well named GB X1. The purpose of this project was the utilization of high temperature tight sedimentary formations for geothermal energy production of about 2 MW thermal to fulfill the energy requirement of surrounding residents [41]. Massive single hydraulic fracturing operation was performed on the targeted formation of Buntsandstein at a depth of 3660 m, by injecting 20,000 m 3 of fresh water carrying no proppants to acquire the fracture area of 1.1 km 2 . The whole fracturing operation was performed in around 5 days with several injection pauses. Several tests were also performed to check stress levels and fracture development and found that the fracture had retained its high conductivity without any proppant usage [41,42].
In our case study, a single hydraulic fracture was used for the verification of the simulated model and then multiple fracturing was performed to analyze multiple fracture response in the presence of stress shadow.

Numerical Simulation Model of Well GB X1
The true representation of the reservoir model was depicted by the actual prevailing geological and strati-graphical conditions. Considering symmetrical geological state, a 1 4 3D model was generated which lies at depth between −3287 and −4100 m and discretized into 36,520 rectangular blocks having dimensions of 1300 (x) × 100 (y) × 813 m (z), respectively ( Figure 13). The injection point is located at a depth of −3660 m in volpriehausen-sandstone formation, while temperature of the model varies according to the geothermal gradient of 0.03 • C/m. The mechanical and hydraulic properties of rocks and stresses were adjusted based on the previous work of Hou et al. [43] and Zhou et al. [44] with the initial pore pressure of about 65 MPa.

Verification of Simulation Model
The study results of Tischer et al. [41] were used for our simulated model verification. A total volume of 20,000 m 3 of pure water was injected with a maximum injection rate of 5.4 m 3 /min. Fracturing operation continued for around 110 h of injection. At the end of the injection, a large fracture having half-length of about 1160 m was created (Figure 14 a). The comparison of the simulated bottom hole pressure (BHP) and measured BHP is shown in Figure 14 b, after the conversion of surface treating pressure to bottom hole pressure, encompassing hydrostatic pressure and frictional losses. The simulated results match well with the measured one, which proves the suitability and acceptability of used parameters and stresses to model fracturing operation. Figure 14 c shows the half fracture geometry at the end of the injection period. The maximum width of 1.85 cm is obtained along the formation having the least value of minimum horizontal stress which also corresponds towards verification of the generated model.

Verification of Simulation Model
The study results of Tischer et al. [41] were used for our simulated model verification. A total volume of 20,000 m 3 of pure water was injected with a maximum injection rate of 5.4 m 3 /min. Fracturing operation continued for around 110 h of injection. At the end of the injection, a large fracture having half-length of about 1160 m was created (Figure 14a). The comparison of the simulated bottom hole pressure (BHP) and measured BHP is shown in Figure 14b, after the conversion of surface treating pressure to bottom hole pressure, encompassing hydrostatic pressure and frictional losses. The simulated results match well with the measured one, which proves the suitability and acceptability of used parameters and stresses to model fracturing operation. Figure 14c shows the half fracture geometry at the end of the injection period. The maximum width of 1.85 cm is obtained along the formation having the least value of minimum horizontal stress which also corresponds towards verification of the generated model.

Results and Discussion of Multiple Fracturing Operation at Well GB X1
After the validation of our geometric/numerical model, an innovative idea is presented to analyze the large geothermal area by creating 10-12 hydraulic fractures to both sides of the well based on the schematic fracture pattern (Figure 15), having a combination of multistage hydraulic fracturing with three vertical wells and lesser injection fluid volumes. This approach can not only provide extra advantages of exploiting a larger geothermal area with lesser seismic activity due to lesser injection volumes, but also reduces the drilling investment cost in comparison with energy production. The accuracy and time consumption for the numerical modeling of multiple fracturing operation is highly dependent on the geometry and number of elements of the simulated model. After conducting numerous simulations, three massive fracture modeling was adopted because of multiple-fracture geometry constraints.

Results and Discussion of Multiple Fracturing Operation at Well GB X1
After the validation of our geometric/numerical model, an innovative idea is presented to analyze the large geothermal area by creating 10-12 hydraulic fractures to both sides of the well based on the schematic fracture pattern (Figure 15), having a combination of multistage hydraulic fracturing with three vertical wells and lesser injection fluid volumes. This approach can not only provide extra advantages of exploiting a larger geothermal area with lesser seismic activity due to lesser injection volumes, but also reduces the drilling investment cost in comparison with energy production. The accuracy and time consumption for the numerical modeling of multiple fracturing operation is highly dependent on the geometry and number of elements of the simulated model. After conducting numerous simulations, three massive fracture modeling was adopted because of multiple-fracture geometry constraints.
The 3D 1 2 geometric model with three injection points is shown in Figure 16a. The horizontal length of the model was reduced to 700 m, while 60 m fracture spacing in y-direction was used as initial base case. The lateral length of the model in the y-direction is 183 m. Through sequential fracturing technique, three fractures were generated after injecting 21,600 m 3 of total fluid volume. Table 2 provides the information about injection parameters and model properties. All three fractures were created by following same injection parameters and conditions with pure water having no proppant. It is observed that minimum horizontal stress has drastically increased from 69 to 95 MPa, especially at the middle of the geometric model (Figure 16b,c).
Energies 2020, 13, x FOR PEER REVIEW 16 of 24 Figure 15. Schematic of the proposed multiple fracture pattern through multi-wells of water injection and production.
The 3D ½ geometric model with three injection points is shown in Figure 16 a. The horizontal length of the model was reduced to 700 m, while 60 m fracture spacing in y-direction was used as initial base case. The lateral length of the model in the y-direction is 183 m. Through sequential fracturing technique, three fractures were generated after injecting 21,600 m 3 of total fluid volume. Table 2 provides the information about injection parameters and model properties. All three fractures were created by following same injection parameters and conditions with pure water having no proppant. It is observed that minimum horizontal stress has drastically increased from 69 to 95 MPa, especially at the middle of the geometric model (Figure 16 b,c).   The 3D ½ geometric model with three injection points is shown in Figure 16 a. The horizontal length of the model was reduced to 700 m, while 60 m fracture spacing in y-direction was used as initial base case. The lateral length of the model in the y-direction is 183 m. Through sequential fracturing technique, three fractures were generated after injecting 21,600 m 3 of total fluid volume. Table 2 provides the information about injection parameters and model properties. All three fractures were created by following same injection parameters and conditions with pure water having no proppant. It is observed that minimum horizontal stress has drastically increased from 69 to 95 MPa, especially at the middle of the geometric model (Figure 16 b,c).   The minimum horizontal stress at the cross section along the x-z plane of each generated fracture provides valuable information related to individual fracture geometry as shown in Figure 17((b)-1,(c)-1,(d)-1). The highest increase in horizontal stress is observed in fracture 3, which is created at the end. The variations in fracture width and shape can be more prominent in massive fracturing operations as higher fluid volumes are used to create large fractures. This can be verified by analyzing the obtained results of fracture geometry with 60 m spacing as illustrated in Figure 17((b)-2,(c)-2,(d)-2). Fracture 1 shows a normal shape having the highest fracture width at the upper section of the formation against the least minimum horizontal stress. The shape, width, and area of 2nd and 3rd consecutive fractures are significantly dependent on the stress re-orientation of 1st and 2nd fractures, respectively. Higher fracture width at the upper side of 1st fracture restricts the propagation of the 2nd fracture at the upper side and forces the propagation of the 2nd fracture at the lower portion. Therefore, maximum fracture width of the 2nd fracture exists at the lower side. In addition, the height of 2nd fracture is increased, while fracture half-length is reduced. It is important to recognize that at greater depth the minimum horizontal stress value is greater as compared to the lower depth (upper side of model), at initial conditions. The maximum width of the 2nd fracture at lower portion is due to the stress re-orientation effect of the 1st fracture, which provides hindrance in the fracture propagation at the upper side.
A similar sort of trend is observed for the 3rd fracture but in the opposite manner. The stress reorientation and fracture width of the 2nd fracture forces the propagation of the 3rd fracture with In addition, the height of 2nd fracture is increased, while fracture half-length is reduced. It is important to recognize that at greater depth the minimum horizontal stress value is greater as compared to the lower depth (upper side of model), at initial conditions. The maximum width of the 2nd fracture at lower portion is due to the stress re-orientation effect of the 1st fracture, which provides hindrance in the fracture propagation at the upper side.
A similar sort of trend is observed for the 3rd fracture but in the opposite manner. The stress re-orientation and fracture width of the 2nd fracture forces the propagation of the 3rd fracture with maximum fracture width at the upper side. Furthermore, configuration of the 3rd fracture is highly distorted with increased fracture half-length as compared to the other two fractures.
Stress shadow effects in the geometric model can be further elaborated by mapping iso-surface contours. Figure 18 demonstrates the comparison of stress shadow iso-surface contours between 60 and 140 m fracture spacing with front, top, and side views. The horizontal well and clip location are highlighted as well. In the case of 60 m spacing, the stress shadow is compacted due to large interference of three fracturing operations, which has created a large distortion in respective fracture shapes. While, in the case of 140 m spacing, contours are fairly distributed in the model representing less intrusion of stress shadow among fracture planes.  To comprehend the stress shadow influence with reference to fracture spacing in this particular area, fracturing simulations were performed further with 80, 120, 140, and 200 m fracture spacing, successively. The corresponding results of three fracture widths along with configurations are shown in Figure 19. To comprehend the stress shadow influence with reference to fracture spacing in this particular area, fracturing simulations were performed further with 80, 120, 140, and 200 m fracture spacing, successively. The corresponding results of three fracture widths along with configurations are shown in Figure 19. The results of 80 m spacing are somewhat similar to 60 m spacing (base case), but with a smooth increase in the 3rd fracture height. Moreover, a maximum width of the 3rd fracture reduces and occurs at the small middle portion of the fracture. This will allow the fracture to propagate in the lower formations as well. As the fracture spacing increases from 80 to 200 m, the homogeneity in the 2nd and 3rd fractures increases. The comparison also shows that the propagation of the 3rd fracture is highly dependent on the configuration of the 2nd fracture. The geometry of the 3rd fracture is somewhat inverse to the 2nd fracture. Therefore, it can be concluded that the continuation of the 2nd and 3rd fracture pattern will occur alternately in the case of a larger number of fracturing operations. The results of 80 m spacing are somewhat similar to 60 m spacing (base case), but with a smooth increase in the 3rd fracture height. Moreover, a maximum width of the 3rd fracture reduces and occurs at the small middle portion of the fracture. This will allow the fracture to propagate in the lower formations as well. As the fracture spacing increases from 80 to 200 m, the homogeneity in the 2nd and 3rd fractures increases. The comparison also shows that the propagation of the 3rd fracture is highly dependent on the configuration of the 2nd fracture. The geometry of the 3rd fracture is somewhat inverse to the 2nd fracture. Therefore, it can be concluded that the continuation of the 2nd and 3rd fracture pattern will occur alternately in the case of a larger number of fracturing operations. Although, with more fracture spacing, chances of similar fracturing patterns increase. However, increased fracture spacing requires a longer horizontal section which can cause an extra burden in drilling, completion, and well operations cost [45]. Considering 140 m spacing as a suitable fracture pattern, the projected 12 fractures pattern is shown in Figure 20. Although, with more fracture spacing, chances of similar fracturing patterns increase. However, increased fracture spacing requires a longer horizontal section which can cause an extra burden in drilling, completion, and well operations cost [45]. Considering 140 m spacing as a suitable fracture pattern, the projected 12 fractures pattern is shown in Figure 20. The summary of this work highlights that appropriate multiple fracture configuration with reference to fracture spacing while incorporating stress shadow effects is essential for assessment of geothermal energy exploitation. Based on this research, geothermal utilization of the North German basin will be conducted in our future work.

Conclusions
In this work, an innovative EGS concept by combining artificially created multiple hydraulic fractures with two horizontal wells for harnessing geothermal energy based on concrete physical performance of the fracture system was presented. The numerical study was conducted using a coupled three-dimensional THM simulator FLAC3D plus and TOUGH2MP-TMVOC to examine the impact of stress shadow on individual fracture fluid flow. Based on this work, the following conclusions are achieved: 1. Through fictive model study outcomes, it is evident that stress shadow superposition affects the subsequent fracture width, which ultimately plays a vital role in geothermal energy production.
The assumption of similar fracture width may lead to erroneous results. 2. Stress shadow superposition enlarges in massive multiple fracturing operations, which eventually distort the fracture propagation as depicted from the results of well GB X1 in the North German basin. 3. Regardless of initial minimum horizontal stress conditions, the shape of the newly created fracture is highly dependent on the previous fracture configuration and attenuation in successive fracture shape becomes more prominent in the case of massive fracturing operations with lower fracture spacing. 4. In order to predict the precise estimations of geothermal energy exploitation, proper multiple fracture configuration becomes imperative and should be studied with reference to fracture spacing under the influence of stress shadow impact. The summary of this work highlights that appropriate multiple fracture configuration with reference to fracture spacing while incorporating stress shadow effects is essential for assessment of geothermal energy exploitation. Based on this research, geothermal utilization of the North German basin will be conducted in our future work.

Conclusions
In this work, an innovative EGS concept by combining artificially created multiple hydraulic fractures with two horizontal wells for harnessing geothermal energy based on concrete physical performance of the fracture system was presented. The numerical study was conducted using a coupled three-dimensional THM simulator FLAC3D plus and TOUGH2MP-TMVOC to examine the impact of stress shadow on individual fracture fluid flow. Based on this work, the following conclusions are achieved: 1.
Through fictive model study outcomes, it is evident that stress shadow superposition affects the subsequent fracture width, which ultimately plays a vital role in geothermal energy production.
The assumption of similar fracture width may lead to erroneous results.

2.
Stress shadow superposition enlarges in massive multiple fracturing operations, which eventually distort the fracture propagation as depicted from the results of well GB X1 in the North German basin.

3.
Regardless of initial minimum horizontal stress conditions, the shape of the newly created fracture is highly dependent on the previous fracture configuration and attenuation in successive fracture shape becomes more prominent in the case of massive fracturing operations with lower fracture spacing.

4.
In order to predict the precise estimations of geothermal energy exploitation, proper multiple fracture configuration becomes imperative and should be studied with reference to fracture spacing under the influence of stress shadow impact. Funding: This research received no external funding.