Geothermal assessment of the Pisa plain, Italy: Coupled thermal and hydraulic modeling

Accepted Manuscript Geothermal assessment of the Pisa plain, Italy: Coupled thermal and hydraulic modeling Guanhong Feng, Tianfu Xu, Fabrizio Gherardi, Zhenjiao Jiang, Stefano Bellani PII:

S0960-1481(17)30341-5

DOI:

10.1016/j.renene.2017.04.034

Reference:

RENE 8729

To appear in:

Renewable Energy

Received Date: 12 December 2016 Revised Date:

10 April 2017

Accepted Date: 15 April 2017

Please cite this article as: Feng G, Xu T, Gherardi F, Jiang Z, Bellani S, Geothermal assessment of the Pisa plain, Italy: Coupled thermal and hydraulic modeling, Renewable Energy (2017), doi: 10.1016/ j.renene.2017.04.034. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT 1 2 3 4

Geothermal assessment of the Pisa plain, Italy: coupled thermal and hydraulic modeling Guanhong Feng a, Tianfu Xu a, Fabrizio Gherardi b,*, Zhenjiao Jiang a, Stefano Bellani b a

5 6

Key Laboratory of Groundwater Resources and Environment, Ministry of Education, Jilin University, Changchun, 130021, China

7 8

b

RI PT

9 10

Istituto di Geoscienze e Georisorse (IGG) – Consiglio Nazionale delle Ricerche (CNR), 56124 Pisa, Italy * Corresponding author. Istituto di Geoscienze e Georisorse, Consiglio Nazionale delle Ricerche, Pisa, Italy. Tel: +39 050 6212316; fax: +39 0503152323. E-mail address: [email protected] (F. Gherardi)

SC

11 Abstract

13

This paper explores the possibility of a development project with a geothermal well doublet in

14

the Pisa plain, Italy. The performance of the system has been evaluated with a 3-dimensional

15

field-scale numerical model that simulates the evolution of temperature and pressure

16

conditions in the aquifer, under different exploitation scenarios. Coupled groundwater flow

17

and thermal transport processes in the reservoir are considered together with non-Darcy fluid

18

flow in the wellbores, and heat exchange between boreholes and surrounding rock formations.

19

Calculations are performed with a parallelized version of the wellbore-reservoir simulator

20

T2Well. This code allows for the efficient modeling of coupled hydraulic-thermal processes

21

over a domain about 40 km2 wide and 1.5 km thick. Simulation results indicate that the energy

22

of the reservoir is sufficient for the designed extraction rate (between 80 and 150 m3/h), but

23

also suitable for much larger rates, up to 250 m3/h. Although aimed at assessing the long-term

24

performance of a specific system, this modeling approach could be profitably applied for the

25

design of similar projects elsewhere.

AC C

EP

TE D

M AN U

12

26 27

Keywords: Low-enthalpy geothermal system; Geothermal doublet; Numerical simulation;

28

Well flow; Aquifer heterogeneity; Pisa plain

29 1

ACCEPTED MANUSCRIPT 30

1. Introduction The European Community has recently promoted initiatives (e.g. [1]) to demonstrate the

32

feasibility of improving local economy and the quality of citizens life through investments in

33

energy efficiency and reduction in carbon emissions. Systemic approaches and organizational

34

innovation are required to achieve by the year 2020 the ambitious target of a 20% reduction of

35

greenhouse gas emissions (compared to 1990 levels) through sustainable production and use of

36

energy [2]. This challenging task calls for an optimized use of different technologies and

37

renewable energy sources, because a single technology/renewable energy could never meet

38

this demand alone. In this context, geothermal energy represents a renewable source with a

39

large potential of energy saving for heating and cooling of buildings.

M AN U

SC

RI PT

31

A sustainable energy action plan that considers the use of geothermal resources for

41

heating/cooling of individual and/or tertiary complexes of buildings has been launched as a

42

part of a smart city pilot project in Pisa, Italy. In this framework, this paper addresses the

43

problem of assessing the feasibility of a geothermal doublet in an area of the Pisa plain (Fig.1)

44

with potential for new urban settlements, and already characterized by the presence of public

45

and private buildings. In this area, refurbishment of existing infrastructures, and planning of

46

new ones, can be easily accomplished. The project relies on a deep geothermal source

47

exploitable at a depth in excess of 600 m b.g.l., and does not consider thermal storage options.

48

Although the energy demand for heating, hot water, air conditioning and cooling cannot be

49

precisely quantified yet, at this stage the project reckons the integration with other renewable

50

sources and the use of electric pumps to provide heating or cooling inside buildings.

AC C

EP

TE D

40

51

The feasibility and the sustainability of a geothermal doublet depend on its power, i.e. on

52

the water yield and the temperature difference between produced and reinjected water, and on

53

its operational lifetime. An economically viable operation requires that a geothermal doublet

54

can be used for a sufficiently long time that depends on site-specific constraints, such as

55

temperature, hydraulic properties of rocks, and the development scheme. Although the

56

theoretical potential of a geothermal power generation is generally assessed under the

57

assumption of a 30 years production life span (e.g. [3,4]), here, the sustainability of a 2

ACCEPTED MANUSCRIPT geothermal doublet has been explored over a larger time span of 50 years. The use of heat

59

energy from medium-temperature geothermal systems for space heating and domestic hot

60

water supply, and the performance of geothermal doublet systems have been extensively

61

investigated since mid-70s’ (e.g. [5-10]). The doublet configuration presents a number of major

62

advantages [11]: (i) it allows for safe reinjection of the produced water; (ii) it improves the

63

heat-recovery factor; (iii) it counteracts production-induced subsidence. Nonetheless,

64

reinjection in geothermal well doublet is expected to develop a thermal front of cold water

65

within the reservoir. Moreover, when there is a contrast in salinity and/or chemical facies

66

between the pristine reservoir water and the reinjected water, a compositional front may be

67

formed as well [12]. Such fronts may be of concern for the sustainability of the geothermal well

68

doublet. In particular, if cold water reaches the production well flowing along the pressure

69

gradient, thermal breakthrough occurs at the producer well, and the temperature of the water

70

extracted progressively declines. Once a certain minimum threshold is reached, the lifetime of

71

the system is sharply shortened. The thermal-breakthrough time depends on geological and

72

hydraulic properties of the reservoir, and on the operational specifics of the doublet.

TE D

M AN U

SC

RI PT

58

In order to evaluate the performance of the Pisa plain thermal system, we used a 3-D

74

numerical model that accounts for the coupled thermal-hydraulic evolution of the local, deep

75

carbonate aquifer exploited with a single geothermal doublet (Fig. 2). The selected numerical

76

simulator allows for solving the challenging problem of the transient non-isothermal,

77

multi-phase, multi-component flow in the integrated wellbore-reservoir system, a problem that

78

could not be efficiently solved by modeling the wellbore or the reservoir separately. Our model

79

also explicitly accounts for the heterogeneous spatial distribution of reservoir properties by

80

exploring a lognormal distribution of permeability, under the assumption of a horizontal

81

aquifer of uniform thickness.

AC C

EP

73

82

At this stage, the assessment of the geothermal potential of the Pisa plain is inherently

83

preliminary, because the size and the type of the district heating/cooling network cannot be

84

defined yet. The main goal of this work is then to build a numerical reservoir model based on

85

key subsurface parameters of the Pisa plain, and to apply it to better understand the 3

ACCEPTED MANUSCRIPT 86

performance sensitivity of a single geothermal doublet. Numerical models provide technical

87

indicators usable in the following phases of the smart city energy action roadmap.

88 2. Site description

RI PT

89

Rock formations in the Pisa plain mainly consist of clastic sediments deposited in a

91

NW-SE striking graben. The area is surrounded by the Pisan Mounts to the north-east, the

92

Livornesi Mounts to the south-west and the Tyrrhenian Sea to the west (Fig. 1). Springs located

93

between the Pisan Mounts and the piedmont plain have temperatures in the range 20-40oC, and

94

a potentially exploitable positive thermal anomaly is observed in the carbonate formation at a

95

depth of 1000-1200 m below ground level [13]. Bellani and Gherardi [14] analyzed the

96

influence of boundary conditions on heat transfer, and outlined the most suitable areas for the

97

geothermal development in the basin. The target layer is represented by a carbonate reservoir

98

more than 1500 m thick. This reservoir is overlain by a low-permeability caprock layer of

99

about 550 m. Two normal faults bound the reservoir on its western and eastern sides.

M AN U

SC

90

The Pisa plain was investigated since the 1970s’ with seismic and gravity surveys for

101

hydrocarbon and geothermal prospection [15]. A few deep (700-3000 m) hydrocarbon

102

exploratory wells were successively drilled in the western part of the basin, reaching a small

103

natural gas reservoir [13], still under exploitation today. Thermal data from deep drillings,

104

along with temperatures from gradient and water wells widely distributed in the area, give

105

geothermal gradients in the range 50-60 °C km-1. The heat flow map shows an average value

106

for the Pisa plain around 100 mW m-2 [14]. In the late 1990s’, a project for the

107

heating/conditioning by geothermal waters of the National Research Council (CNR) campus in

108

Pisa was set up. The project foresaw the drilling of two deviated wells, in the 800-1200 m depth

109

range, for fluid withdrawal and reinjection. Only one well (“S.Cataldo 1” well) was drilled at

110

that time, down to about 1050 m (850 m vertical). A temperature of 49.5 °C was measured in

111

water coming from the top of the carbonate reservoir, at a depth of about 600 m, during a

112

pumping test. Chemical analyses of major constituents and in situ measurements of

113

physico-chemical parameters on water samples, revealed near-neutral pH, medium salinity (up

AC C

EP

TE D

100

4

ACCEPTED MANUSCRIPT to 10 g/kgw), a Na-Cl-HCO3 facies, and negligible gas contents [16]. This suggested that

115

formation waters are not chemically aggressive, making the area suitable for geothermal

116

development. Due to unfavorable economic conditions and technical difficulties to fit the

117

geothermal resource to an already installed natural gas fired plant, the well was then abandoned,

118

giving up the whole project. The temperature and delivery data though showed a good potential

119

for further geothermal investigations in the area.

RI PT

114

120

Figure 1. Location map of the area under study in the Pisa plain, with indication of outcropping geological formations and their permeability. The blue contour is the approximate extension of the

SC

121 122 123

modelled area.

125

M AN U

124

Figure 2. Setup of the thermal system and geological framework.

126 127

3. Model setup

Laterally constrained by the faults, the modelled domain has a length of 8000 m on the top,

129

and 9050 m on the bottom, along the east-west direction (Figs. 3a and 3b), significantly larger

130

than along the north-south direction (5000 m). The domain has a vertical extent of 1500 m

131

downwards from the top of the carbonate reservoir, large enough to eliminate the influence of

132

boundary conditions at the bottom of the system.

EP

TE D

128

Following previous studies [17], our reservoir model considers a specific doublet

134

configuration, that exploits two directional wells, drilled from the same drill pad, with a

135

distance of 1100 m between well bottoms at a depth of 1000 m b.g.l. (Fig. 3c). The wells are

136

symmetrically inclined, and located in the central part of the modelled volume (Fig. 3a). In the

137

upper 200m, both wells are vertically drilled through the caprock, with a diameter of 0.7 m.

138

The diameter decreases to 0.5 m between 200 to 550m, where an inclination of 15° is

139

assumed. Within the carbonate formation, the inclination further increases with an angle of

140

45o between the wellbore and the vertical. Both boreholes are slotted from 550 to 1000 m

141

b.g.l. in correspondence to the carbonate formation, where they penetrate for 450 m with a

142

diameter of 0.3 m. The total length is 1200 m, equivalent to a vertical penetration depth of

AC C

133

5

ACCEPTED MANUSCRIPT 1000 m. In this configuration, the positioning, geometry and spacing of the wells respond to a

144

trade-off between efficiency and operational lifetime of the doublet. In fact, wells should not be

145

placed too far apart to maximize the main benefit from reinjection (i.e. pressure sustain), and

146

not too close to avoid thermal short-circuiting. The optimal positioning and spacing of a

147

geothermal doublet is largely dependent on the economic lifetime of the project, and

148

configurations different from that explored here could be taken into consideration in future

149

developments, possibly according to updated societal demands and economic constraints.

RI PT

143

150

Figure 3. Lateral 2D cross-section of the domain (a) and 3D view of the model grid (b), together with

SC

151 152

wellbore structure of the geothermal doublet (c).

154

M AN U

153 3.1. Numerical model

The coupled thermo-hydraulic processes in the reservoir, and the aqueous flow in the

156

wellbores are modelled with the fully coupled wellbore-reservoir simulator T2well [18].

157

T2well was developed by introducing wellbore flow equations into the reservoir simulator

158

TOUGH2 [19]. In this code, wellbores and reservoir are handled as separate sub-domains.

159

The flow process in the reservoir is mathematically described by mass balance equations

160

based on Darcy's law, while the flow in the wellbores is governed by 1-D momentum

161

conservation equation [20,21]. The equations for mass and energy balance are given in

162

Appendix A. The water exchange between wellbore and reservoir is expressed by Darcy's law.

163

Based on the momentum equation, the drift-flux model is applied to calculate the velocity of

164

each fluid phase. For this study, a parallel version of T2well (T2well-MP) has been used to

165

overcome the problem of the huge computational burden associated with the large number of

166

grid cells considered (more than 105 grid cells in total), and the coupled hydraulic and

167

thermal processes intervening in the numerical model at the well-reservoir scale. T2well has

168

been parallelized by following the same procedures used for the parallelization of the

169

TOUGH2 code (TOUGH-MP code [22,23]). Efficient calculations are guaranteed by a

170

specific task allocation for the processors that avoids the exchange of data among wellbore

171

cells and minimizes the exchange between wellbore and reservoir grid cells. Full details on

AC C

EP

TE D

155

6

ACCEPTED MANUSCRIPT 172

equations and flow chart of the parallelized T2well simulator are given in Appendix B.

173 174

3.2. Initial and boundary conditions The model domain is laterally constrained by two NW-SE trending low-permeability

176

faults [14], and is overlain by a 550 m thick succession of clay and sand mixed layers that

177

acts as an impermeable caprock. Therefore, the lateral and the top surfaces are defined as

178

isolated boundaries. The aquifer is horizontal, mineralogically homogeneous (carbonate

179

rocks), and of uniform, considerable thickness (>2000m) over the area of interest. The

180

thickness of the simulated system is 1500 m, which implies that the bottom of the domain is

181

sufficiently far away from the doublet to approach a condition of no water flow boundary

182

(Fig. 3). Heat flow is allowed from this boundary, and it is dynamically calculated with a

183

semi-analytical equation [24], in agreement with the temperature gradient measured below

184

the depth of 575 m. The setup of different boundary conditions (e.g. constant head and no

185

flow boundary conditions at fault positions) has a negligible effect on the numerical outputs.

M AN U

SC

RI PT

175

The temperature measured at the depth of 575 m (the top of carbonate reservoir) is 49oC

187

[14], whereas the temperature of the shallow aquifer is assigned as 15°C, the mean local

188

ambient temperature. The caprock is characterized by a vertical gradient of about 59oC/km.

189

Within the reservoir, the temperature gradient falls to about 15oC/km. Density, porosity and

190

thermal conductivity of carbonate reservoir rocks are taken from the literature (e.g. [25,26])

191

and summarized in Table 1.

193 194

EP

AC C

192

TE D

186

Table 1. Thermal and hydrological parameters of wellbores and carbonate reservoir

195

The permeability is assigned according to the results of an injection test performed in

196

1998 in an exploratory well located inside the present study area (San Cataldo 1 well;

197

CNR-IGG unpublished data). During this test, the static pressure at wellhead kept 0.2MPa.

198

The injectivity was estimated at 250 to 300 m3/hour/MPa. The pumping rate maintained 108 7

ACCEPTED MANUSCRIPT m3/hour during the 33 hours of the test. As a result, the models consider a carbonate reservoir

200

permeability of 5.0 × 10-14 m2, and a vertical permeability of 5.0 × 10-15 m2 is assumed with

201

an arbitrary horizontal to vertical permeability ratio of 10.

202

Under the conditions above, we compute the temperature evolution in the reservoir over a

203

50-years operational lifetime, by envisaging two possible scenarios, both energetically

204

plausible considering the mild climate of the area (15°C annual average air temperature): (a)

205

a constant pumping rate (80 m3/hour) throughout the year, and (b) a higher pumping rate (200

206

m3/hour) limited to the winter months (November – March) and no production for the

207

remaining seven months of the year. The sensitivity to different injection/extraction rates is

208

also explored over an 80 to 600 m3/hour interval with the aim to define a maximum

209

sustainable production rate. All these scenarios consider no surplus water (extraction rate =

210

injection rate). Moreover, the fluid injection temperature is assumed at 15°C, the local

211

average groundwater temperature. This is an approximated reference value, but different

212

reinjection temperatures are possible, depending on the adopted operational scheme.

M AN U

SC

RI PT

199

TE D

213 4. Results and discussion

215

4.1. Homogeneous reservoir

216

4.1.1. Constant rate

EP

214

A reference case simulation (RCS), with a constant injection/production rate of 80 m3/hr

218

over 12 months per year, is run by assuming homogeneous hydrogeological and thermal

219

parameters (Fig. 4). Over the long-term, this configuration induces minor variations in the

220

hydraulic pressure in all the reservoir, included the sector near the production well (Fig. 4b).

221

In this part of the domain, the largest pressure drop computed is in fact in the order of 1 MPa

222

or lower. Conversely, pressure slightly increases around the production well. Overall,

223

groundwater dynamics is controlled by the downward movement of cold water near the

224

injection well and by the upward movement of warm water near the production well. In

225

particular, the steady injection of 15°C water causes temperature to go below 30°C in a large

226

volume of rocks, down to about 500 m below the top of the reservoir. Noteworthy, the

AC C

217

8

ACCEPTED MANUSCRIPT migration of cold injected water does not lower the outflow temperature of the water

228

delivered by the extraction well, testifying for a suitable design of the geothermal doublet,

229

under the conditions considered in the homogeneous reservoir model. This is because the

230

downward migration of cold water does not result in any short-circuiting between producer

231

and injector, but causes the onset of a forced convection that drives more thermal water from

232

deeper sectors of the reservoir towards the extraction well. As a result, after 50 years, the

233

temperature increases by about 2.4±1.2°C in the production well (bottom: +1.2°C; top:

234

+3.6°C).

RI PT

227

Figure 4. Temperature (a) and pressure (b) distribution after 50 years along a 2D WSW-ENE trending section of the modelled domain (reference case simulation, RCS, with 80 m3/h injection/production constant rate).

239

M AN U

236 237 238

SC

235

To explore the geothermal potential of the reservoir, we investigate some additional cases

241

with both the pumping and the injection rates increasing through several steps from 80 to 600

242

m3/hour. For instance, under the 600 m3/hour injection/production rate condition (Fig. 5), a

243

quite large cold plume (T<25°C) is predicted to develop into the reservoir, around the

244

production well, over a volume having an approximate diameter of 1.5 km, and a vertical

245

extent of about 0.6 km. Injection/extraction rates higher than 150 m3/hour could result in

246

measurable wellhead temperature drops, whereas heat extraction rates do not show any

247

significant decrease (over 50 years of operation) up to about 250 m3/hour. For example, under

248

the 400 and the 600 m3/hour conditions, the temperature starts decreasing after less than 10

249

years, with an overall predicted reduction of 6 and 8°C after 50 years, respectively (Fig. 6a).

250

Under the same conditions, the heat extraction rate is predicted to decrease by 10% (from 15.5

251

to 14 MW at 400 m3/hour) and 24% (from 25 to 19 MW at 600 m3/hour), respectively. The 150

252

m3/hour condition is then cautiously considered as an operational limit along a hypothetically

253

maximum 50 years lifespan of the geothermal doublet, under the homogeneous reservoir

254

assumption.

AC C

EP

TE D

240

255 9

ACCEPTED MANUSCRIPT 256 257

Figure 5. Temperature distribution after 50 years along a 2D WSW-ENE trending section of the modelled domain (400 m3/hour case). Red dashed lines are the traces of the wells.

258 The patterns for heat extraction rate (Fig. 6b) mimic the temperature behavior because the

260

mass flow is constant during the calculation, and the pressure has a negligible effect on energy

261

balance. The heat extraction rate (G) is calculated in fact with the following equation:

262

(1)

RI PT

259


= () × ℎ() − ( ) × ℎ( )

265

() and () stand for the production and injection well, respectively. The different

operational regimes investigated in our simulations are predicted to extract different amounts

266

of heat from the rocks, with production plateaus quickly achieved in the early stages of

267

operation (i.e. few months to few years) in the range 80 to 400 m3/hour. The highest the

268

extraction/injection rates (>400 m3/hour), the largest the decrease in wellhead temperature and

269

heat extraction rate, and the induced pressure drawdown along the production well (up to -0.45

270

MPa in a radius of about 100 m from the well axis; Fig. 7).

275 276 277 278

M AN U

TE D

274

Figure 6. Temperature (a) and heat extraction rate (b) variation for different constant production rates.

Figure 7. Pressure field after 0.1 and 50 yrs (reference case simulation, RCS, with 80 m3/h

EP

272 273

injection/production constant rate).

AC C

271

SC

264

where  is the mass flow rate [kg s-1], ℎ the specific enthalpy [kJ kg-1], and the subscripts

263

4.1.2. Intermittent production-recovery cycles

279

This scenario considers injection/production to occur over a period of five months per

280

year at a constant rate of 200 m3/hour, with stop of production during the remaining seven

281

months of the year. The total thermal water outflow (around 720 × 103 m3/year) is almost the

282

same as the volume extracted in the RCS (constant rate of 80 m3/hour over 12 months per

283

year). The temperature and pressure distribution in the reservoir is similar to the RCS.

284

However, under the 5-month operation regime (Scenario 2), the wellhead temperature is 10

ACCEPTED MANUSCRIPT 285

0.5°C higher than for the RCS (Scenario 1; Fig. 8). Noteworthy, the model predicts a quick

286

thermal recovery, with full restoration of the initial aquifer temperature during the 7-month

287

period of no-production.

288 Figure 8. Time variation of the outflow temperature for two different scenarios. “Scenario 1”: 80 m3/hour constant rate over 12 months per year; “Scenario 2”: 200 m3/hour constant rate over 5 months per year. “Scenario 2” curves plot only data for the 5 month period of production. During the remaining 7 months without production (not shown), temperature goes down to 20°C, and quickly rises to more than 50°C as soon as production re-starts.

RI PT

289 290 291 292 293

295

SC

294 4.1.3. Effect of wellbore-reservoir coupling

The effect of wellbore-reservoir coupling is evaluated with a simulation that considers

297

the same initial and boundary conditions of the RCS, but neglects temperature and pressure

298

variations predicted by T2well along the production well (“uncoupled model”; Fig. 9). It can

299

be seen that the uncoupled model overestimates roughly by 10% (about +0.2 to +0.3 MW)

300

the heat extraction rate of the geothermal doublet, because it neglects temperature (heat loss

301

along the pipe) and pressure (frictional pressure drops) variations expected within the

302

wellbore. These variations are accounted only by the coupled model. The wellbore-reservoir

303

coupled configuration then emerges as a necessary approach to evaluate the feasibility of the

304

geothermal doublet and to increase the reliability of its performance assessment.

308 309 310

TE D

EP

306 307

Figure 9. (a) Heat extraction rates predicted by the coupled and uncoupled models in the RCS . (b) Temperature profile along injection and production wells after 50 years (coupled model).

AC C

305

M AN U

296

4.2. Heterogeneous reservoir

311

The heat extraction from rocks, and the thermal recovery, depend on the structure of the

312

flow-paths in the reservoir. Natural aquifers are intrinsically heterogeneous, and aquifer

313

heterogeneities are expected to influence fluid flow, thermal retardation, and then the

314

profitability of geothermal well doublets. In the case of Pisa plain, the impact of 11

ACCEPTED MANUSCRIPT 315

heterogeneities on the thermal behavior of the reservoir is investigated with the aim to identify

316

the conditions for possible thermal-breakthroughs and related temperature decline. In this

317

numerical exercise, permeability acts as a master variable, because the assumption is made that

318

changes in permeability have the largest influence on thermal lifetime and pressure drop. Based on results from calculations with homogeneous distribution of hydrological and

320

thermal parameters, a set of additional models is run by assuming the permeability to be: (i)

321

lognormally distributed [27,28] in a sub domain (heterogeneous volume, hereinafter HV)

322

identified by the following coordinates: X = 6500±1000 m, Y = 2500±500 m, Z = 500±500 m;

323

(ii) homogeneously distributed (Table 1) over the rest of the domain (Fig. 3). In the setup of

324

these models, permeability is considered to behave as a regionalized variable within the HV.

325

Standard geostatistical techniques are then applied to this sub-domain by assuming an uniform

SC

RI PT

319

328

values of 5.0 × 10-12 m2 and 5.0 × 10-16 m2, respectively. Permeability values are finally

329

assigned to the 16000 grid cells of the ECG by calculating a spherical semivariogram (Fig. 8)

333 334

335

336

TE D

332

without nugget effect ( = 0), with a sill of 0.8 ( = 0). The correlation length strongly

depends on the variogram type and the model scale (e.g. between 0.5 to 800m [29-31]). Here, a correlation length of 500 m ( = 500) is chosen. Overall, the following relationships hold:

EP

331

of 5.0 × 10-14 m2, and a variance of 0.8 [29] , as to obtain maximum and minimum permeability

(ℎ) = × #$(%) − (% + ℎ)' = × ( × $(%) − (% + ℎ)'" " " !

AC C

330

M AN U

327

grid spacing (∆ = ∆ = ∆ = 50; ”external cubic grid”, ECG), a mean permeability value

326

(ℎ) = ) +  × +

,-

".

0, ℎ = 0

!

− "./ 0 , 0 < ℎ ≤  -/

 + , ℎ ≥ 

(2)

(3)

337

The semivariogram allows for a quantitative description of the spatial dependence of a

338

randomly heterogeneous field [32,33]. It is defined as the variance of the difference between

339 340

any two spots in the field, as given by equation (2). In this framework, % is a random spot in the field, ℎ is the distance from % to any other spot,  is the spatial variable (permeability in 12

ACCEPTED MANUSCRIPT

342

our case), (ℎ) is the correlation of any two spots, and higher (ℎ) values indicate weaker

343

lognormal distribution of values, and those generated in our model by applying a geostatistical

344

approach is given in Figure 10. The initialization of the model is finalized by passing the

345

permeability values calculated over the ECG to the HV sub-domain of the field-scale model,

346

correlation degrees. A comparison between the permeability values obtained by assuming a

that has a different spatial discretization (∆ = ∆ = 100 ; ∆ = 50 ). This is done by

RI PT

341

347

adjusting (by kriging) the values computed over the ECG to the grid of the field-scale model

348

(HV sub-domain).

SC

349

Figure 10. (a) Lognormal distribution of permeability (variance: 0.8; mean: 5.0 × 10-14 m2. (b) Experimental semivariogram (spherical model) inferred from 100 generated permeability values.

353

The 7 − 8 relation in geologic media depends on the intrinsically complex structure of

354

M AN U

350 351 352

the porous space, and in particular on the dimensions of connected pores/fractures (e.g.

356

[34-36]). Although experimental observations and field logs have demonstrated that no single

357

model can reliably describe this intrinsic complexity, different types of models have been

358

proposed to link permeability and porosity values (e.g. [37-40]). Based on the qualitative

359

observation that permeability values generally correlate with porosity values, here, for the

360

sake of simplicity, we assume that these variables obey the following equation [34]:

EP

361

TE D

355

log(8) =  × 7 + <

(4)

364

average and maximum porosity values calculated with this equation for the Pisa plain

365

carbonate aquifer are 0.01, 0.08 and 0.15, respectively.

AC C

363

where 8 is the permeability (m2) and 7 is the porosity,  and < are two empirical

362

coefficients set to 28.5 and 15.9 in the case of carbonate reservoirs [41,42]. The minimum,

366

Infinite types of heterogeneous permeability fields can be generated by means of random

367

number generators. In the following, only three selected cases are discussed as representative

368

of these limiting conditions: (i) high-permeability zone between the two wells (Figs. 11a and

369

11b; heterogeneous case 1, hereinafter HC1); (ii) low-permeability zone between the two 13

ACCEPTED MANUSCRIPT wells (Figs. 11c and 11d; heterogeneous case 2, HC2); (iii) low-permeability zone at the

371

bottom of the injection well and high-permeability zone at the top of the reservoir (Figs. 11e

372

and 11f; heterogeneous case 3, HC3). Although these are only a few possible realizations of

373

the subsurface structure, and, therefore, they are intrinsically not representative for any part

374

of the reservoir, they are used to explore the sensitivity of the numerical model and to

375

quantify the effects of heterogeneities on the profitability of a geothermal project with the

376

limits in porosity and permeability so far defined.

RI PT

370

Figure 11. Three hypothetical heterogeneous permeability distributions (a: HC1, c: HC2, e: HC3) together with the corresponding 2D cross-sections drawn at Y=2500 (b, d, f). Red dashed lines provide the approximate location of the injection and production wells.

M AN U

378 379 380

SC

377

381

The models used to investigate the effects of heterogeneity consider the same rate of the

383

RCS (constant 80 m3/h, 12 months per year). The temperature profiles along the production

384

wellbore (Fig. 12) highlight the importance of (i) the high-permeability zones at the top of the

385

reservoir (HC3) and between the two wells (HC1) in conveying low temperature fluids

386

towards the production well, and of (ii) the low-permeability zone between the two wells

387

(HC2) in contrasting any tendency of fluid short-circuiting between the injector and the

388

producer.

390 391 392

EP

AC C

389

TE D

382

Figure 12. Temperature profile after 50 years in the production well in the reference case simulation (RCS) and in three selected heterogeneous cases (HC1, HC2, HC3).

393

As expected, the cold plume associated with the reinjection process preferentially

394

propagates within the reservoir in those areas where the distribution of permeability values

395

allows for efficient fluid circulation around the injection well (Fig. 13). Accordingly, the HC1

396

and HC3 models predict a large migration of cold water (t<35°C) along the vertical direction, 14

ACCEPTED MANUSCRIPT down to -800 m below the top of the reservoir, and eastwards, up to about 6700 m from the

398

western boundary of the reservoir. The HC2 model mimics a completely different scenario,

399

where the low permeability zone between the wells acts as a barrier for the migration,

400

allowing the extraction well to efficiently deliver local thermal fluids over the 50 year

401

lifespan of the geothermal doublet. In this case, the effect of the barrier is not only to prevent

402

the short-circuiting of reinjected water, but also to enhance production from deeper and

403

warmer levels of the aquifer. This favorable condition is highlighted by higher temperature

404

values of the delivered fluid, compared to the RCS (Fig. 12). Despite this evidence, there is

405

no direct and univocal correlation between the enhanced productivity of the geothermal

406

doublet and the decreased permeability of the rocks between the injector and the producer.

407

Too low permeability values might cause in fact the reverse, undesired effect of decreasing

408

the injection and production flow rates, due to the reduced efficiency of flow paths between

409

the two wells. The heterogeneous distribution of permeability causes the geothermal doublet

410

to differently perform with depth in terms of extracted energy (Fig. 14a) and production

411

rates/pressure (Fig. 14b). Compared to the RCS, the HC1 and HC3 models predict in fact that

412

most of the thermal supply is provided by layers near the top of the reservoir, and the reverse,

413

the HC2 model predicts production to occur at the bottom of the wellbore, with almost no

414

production in the upper part of the reservoir (Fig. 14b). The HC2 model emerges as the most

415

favorable scenario also in terms of amounts of extracted energy, with an almost constant rate

416

of 3.6 to 3.7 MW over the investigated lifespan of the geothermal doublet (Fig. 14a). The

417

RCS wellhead pressures are the lowest for the production, and the highest for the injection

418

well, respectively (Fig. 15).

AC C

EP

TE D

M AN U

SC

RI PT

397

419

Independently from the permeability field assumed in the calculations, no significant

420

pressure drawdowns are calculated over the 50 years lifespan of the geothermal doublet. This

421

is interpreted as an effect of the elevated transmissivity of the aquifer. Overall, the geothermal

422

doublet of the Pisa plain is predicted to operate in a sustainable manner over decades without

423

significant production temperature drop. Although this conclusion is specific to the hydraulic

424

conditions and the operational scenarios considered by the numerical models, the predicted 15

ACCEPTED MANUSCRIPT scenarios of the production capacity over time can be considered reliable since the

426

simulations are calibrated with results from preliminary hydraulic tests, and the models are

427

initialized with geological and hydrogeological information specific to the site. An intrinsic

428

limitation of the present analysis is that it does not take into account the physical and

429

chemical processes linked to the formation, migration and accumulation of mineral scales

430

and/or deposits in the wells (mechanically transported fine particles) and/or in the reservoir

431

(chemically-induced mineral precipitation). Under particular circumstances, these processes

432

have the potential to clog and damage boreholes, pipelines and source rocks (e.g. [43,44]).

433

The amount and kind of minerals that might precipitate, and/or the amount of solid particles

434

that may accumulate as an effect of the perturbation induced in the aquifer by the geothermal

435

doublet, deserve thus further attention. A geochemical and reactive transport analysis of the

436

processes induced in the aquifer by the geothermal doublet is then needed for a more reliable

437

prediction of the sustainability of the heat production from the carbonate aquifer over several

438

decades. This work has not been done yet due to the lack of detailed information on local

439

mineralogy and chemical composition of pore and reinjected waters.

M AN U

SC

RI PT

425

Figure 13. 2D WSW-ENE section showing the temperature distribution after 50 years in a homogenous reservoir (a: RCS) and in three hypothetical heterogeneous configurations of the same reservoir (b: HC1; c: HC2; d: HC3). Red dashed lines indicate the approximate location of the injection and production wells.

EP

441 442 443 444

TE D

440

445

Figure 14. Heat extraction (a: variation over the 0 to 50 years interval) and production flow (b: snapshot after 50 years) rates in a homogenous reservoir (RCS) and in three hypothetical heterogeneous configurations of the same reservoir (HC1, HC2, HC3).

450 451 452

Figure 15. Time variation of wellhead pressure in the production (a) and injection (b) well. Four

AC C

446 447 448 449

scenarios are compared: homogenous reservoir (RCS) vs. three hypothetical heterogeneous configurations of the same reservoir (HC1; HC2; HC3).

453 454 455

5. Conclusions Because of its ability to ensure reliable and affordable heating/cooling by means of a 16

ACCEPTED MANUSCRIPT low-carbon energy carrier, space heating by geothermal technology is among the proven

457

renewable energy solutions considered in the planning of future smart cities. As a very first

458

step of a wider smart city energy action plan, the performance of a geothermal doublet

459

thermal system in the neighborhood of Pisa, Italy, is evaluated with a parallelized version

460

(T2well-MP; [11-12]) of the T2well simulator [7]. This code allows for the efficient modeling

461

of the transient non-isothermal, multi-phase, multi-component flow in the integrated

462

wellbore-reservoir system, over a large grid of more than 105 grid cells.

RI PT

456

Based on preliminary pumping tests and information of hydraulic and petrophysical

464

properties of local rocks, numerical simulations consider a geothermal doublet that exploits

465

two directional wells penetrating the local carbonate reservoir down to a depth of 1000 m

466

below ground level. The sustainability of the geothermal doublet is explored over a time span

467

of 50 years, accounting also for the effects of three selected heterogeneous permeability fields

468

on thermal regime and production rates. Heterogeneities significantly affect production

469

temperatures and rates, and advantageous doublet configurations are identified by sensitivity

470

analysis. The carbonate reservoir in the Pisa plain turns out to have an elevated transmissivity

471

and to contain sufficient resource for the designed extraction rates of 80 m3/hour over 12

472

months per year (continuous injection/production rate), or 200 m3/hour over 5 months per year

473

(intermittent injection/production rate). No cold front breakthroughs are predicted, testifying

474

for the suitable design of the geothermal doublet.

EP

TE D

M AN U

SC

463

Overall, our numerical models allow for a realistic description of the processes occurring

476

at the injection and production wells, and set the baseline for an assessment of the long term

477

behavior of a geothermal doublet in the Pisa plain. These results should be considered as a

478

preliminary support tool for the location and the dimensioning of a geothermal doublet.

479

Further, detailed information about water chemistry, reservoir mineralogy and petrophysical

480

parameters are needed to run geochemical and reactive transport models aimed at evaluating

481

the long term impact of water quality on geothermal exploitation. In addition, economic

482

factors, not considered in the present work, should be taken into account for an optimized

483

configuration of the geothermal facility, because doublets should be operated at a rate

AC C

475

17

ACCEPTED MANUSCRIPT 484

corresponding to the installed capacity of their heating capabilities. More accurate key

485

parameters data for the chosen technology will be necessary in the next design phase of the

486

district heating system, to size the surface system to achieve optimum performance with the

487

lowest possible cost.

489

RI PT

488 Acknowledgements

This work was jointly supported by the National Natural Science Foundation of China

491

(Grant No. 41572215) and by the 111 project (No. B16020). Two anonymous reviewers are

492

warmly acknowledged for their useful comments.

SC

490

494 495

References

496

[1]

498

Energy

Europe,

Support

for

sustainable

energy

solutions;

http://ec.europa.eu/energy/intelligent (Accessed 7 December 2016) [2]

Energy

Performance

of

Buildings

TE D

497

Intelligent

M AN U

493

Directive

499

http://eur-lex.europa.eu/legal-content/en/TXT/?uri=celex%3A32010L0031

500

December 2016)

2010/31/EU; (Accessed

7

[3] C.F. Williams, M.J. Reed, R.H. Mariner, A review of methods applied by the U.S.

502

Geological Survey in the assessment of identified geothermal resources, USGS Open File

503

Report 1296 (2008) 30 pp.

505

AC C

504

EP

501

[4] G. Beardsmore, L. Rybach, D. Blackwell, C. Baron, A protocol for estimating and mapping the global EGS potential, GRC Transactions 34 (2010) 301-312.

506

[5] A.C. Gringarten, J.P. Sauty, A theoretical study of heat extraction from aquifers with

507

uniform regional flow, J. Geophys. Res. 80 (1975) 4956-4962.

508

[6] C.F. Tsang, M.J. Lippmann, P.A. Whiterspoon, Production and reinjection in geothermal

509

reservoirs. GRC Transactions 1 (1977) 301-303.

18

ACCEPTED MANUSCRIPT [7] C.R. Faust, J.W. Mercer, Geothermal reservoir simulation. 2. Numerical solution

511

techniques for liquid- and vapour-dominated hydrothermal systems, Wat. Resour. Res. 15

512

(1979) 31-45.

513

[8] M.J. Lippmann, C.F. Tsang, Groundwater use for cooling: associated aquifer temperature

514

changes, Ground Water 18 (1980) 452-458.

515

[9] J.W. Mercer, C.R. Faust, W.J. Muller, F.J. Person Jr, Review of simulation technologies

516

for aquifer thermal energy storage (ATES), Adv. Hydrosci. 13 (1982) 1-129.

517

[10] D. Banks, An introduction to thermogeology: ground source heating and cooling,

518

Blackwell, Oxford, UK (2008) 339 pp.

519

[11] K. Bucher, I. Stober, Geothermal energy: from theoretical models to explorations and

524 525 526 527 528 529 530

SC

M AN U

523

for urban heating, Pure Appl. Geophys. 117 (1979) 297-308.

[13] S. Bellani, S. Grassi, P. Squarci, Geothermal characteristics of the Pisa plain, Italy. In: Proc World Geothermal Congress, Florence, Italy, 2 (1995) 1305-1308. [14] S. Bellani, F. Gherardi, Thermal features of the Pisa plain, a Neogenic basin in central

TE D

522

[12] A.C. Gringarten, Reservoir lifetime and heat recovery factor in geothermal aquifers used

Italy, GRC Transactions 38 (2014) 357-361.

[15] M. Mariani, R. Prato, I bacini neogenici costieri del margine tirrenico: approccio sismico-stratigrafico, Mem. Soc. Geol. It. 41 (1988) 519-531 (in Italian).

EP

521

development, Springer, Berlin & Heidelberg (2013).

[16] S.Biagi, F.Gherardi, G.Gianelli, A simulation study of CO2 sequestration in the Arno River Plain (Tuscany, Italy), Energy Sources Part A 28 (2006) 923-932.

AC C

520

RI PT

510

531

[17] S. Bellani, G. Buonasorte, S. Grassi, P. Squarci, Geological and structural features of the

532

San Cataldo district heating project (Pisa, Italy), Proceedings 5th World Renewable Energy

533

Congress, Florence 4 (1998) 2746-2749.

534 535 536 537

[18] L. Pan, C.M. Oldenburg, T2Well - An integrated wellbore–reservoir simulator. Comput. Geosci. 65 (2014) 46-55. [19] K. Pruess, C.M. Oldenburg, G. Moridis, TOUGH2 User's guide, version 2.0. LBNL Report 43134, Berkeley, CA, 1999.

19

ACCEPTED MANUSCRIPT 538

[20] L. Pan, C.M. Oldenburg, Y-S. Wu, K. Pruess. T2Well/ECO2N version 1.0: multiphase

539

and non-isothermal model for coupled wellbore-reservoir flow of carbon dioxide and

540

variable salinity water. LBNL Report 4291E, Berkeley, CA, 2011.

542

[21] L. Pan, S.W. Webb, C.M. Oldenburg, Analytical solution for two-phase flow in a wellbore using the drift-flux model, Adv. Water Resour. 34 (2011) 1656-1665.

RI PT

541

543

[22] K. Zhang, Y-S Wu, G.S. Bodvarsson, Parallel computing simulation of fluid flow in the

544

unsaturated zone of Yucca Mountain, Nevada, J. Contam. Hydrol. 62-63 (2003) 381-399.

545

[23] Zhang K, Y-S. Wu, K. Pruess, User's guide for TOUGH2-MP - a massively parallel

548

SC

547

version of the TOUGH2 code. LBNL Report 315E, Berkeley, CA, 2008.

[24] P. Vinsome, J. Westerveld, A simple method for predicting cap and base rock heat losses in thermal reservoir simulators, J. Can. Petrol. Technol. 19 (1980) 87-90.

M AN U

546

[25] P.A. Domenico, F.W. Schwartz, Physical and chemical hydrogeology, John Wiley and

550

Sons, New York (1990) 824 pp.

551

[26] J.P. Brill, H. Mukherjee, Multiphase flow in wells, Society of Petroleum Engineers, SPE

552

Monograph Series (1999).

553

[27] S. Haldar, G.L.S. Babu, Effect of soil spatial variability on the response of laterally loaded

555 556

pile in undrained clay, Comput. Geotec. 35 (2008) 537-547. [28] G.B. Baecher, J.T. Christian, Reliability and statistics in geotechnical engineering, John Wiley & Sons, London and New York (2003).

EP

554

TE D

549

[29] H. Tian, F. Pan, T. Xu, B.J. McPherson, G. Yue, P. Mandalaparty, Impacts of hydrological

558

heterogeneities on caprock mineral alteration and containment of CO2 in geological storage

559

sites, Int. J. Greenh. Gas Control 24 (2014) 30-42.

AC C

557

560

[30] A. Srivastava, G.L.S. Babu, S. Haldar, Influence of spatial variability of permeability

561

property on steady state seepage flow and slope stability analysis, Eng. Geol. 110 (2010)

562

93-101.

563

[31] A.W. Western, S-L. Zhou, R.B. Grayson, T.A. McMahon, G. Blöschl, D.J. Wilson,

564

Spatial correlation of soil moisture in small catchments and its relationship to dominant

565

spatial hydrological processes, J. Hydrol. 286 (2004) 113-34. 20

ACCEPTED MANUSCRIPT

567 568 569 570 571 572 573

[32] P.J. Curran, The semivariogram in remote sensing: an introduction. Remote Sens. Environ. 24 (1988) 493-507. [33] D.L. Zimmerman, M.B. Zimmerman, A comparison of spatial semivariogram estimators and corresponding ordinary kriging predictors. Technometrics 33 (1991) 77-91. [34] G.E. Archie, Introduction to petrophysics of reservoir rocks, AAPG Bull. 34 (1950)

RI PT

566

943-961.

[35] G.E. Archie, Classification of carbonate reservoir rocks and petrophysical considerations, AAPG Bull. 36 (1952) 278-298.

[36] J. Bear, Dynamics of fluids in porous media, Dover Publications, New York (1972).

575

[37] P.C. Carman, Flow of gases through porous media, Academic Press, New York (1956).

576

[38] A. Timur, An investigation of permeability, porosity, and residual water saturation

579 580 581

M AN U

578

relationships for sandstone reservoirs, The Log Analyst 9 (1968) 8-17. [39] A.J. Katz, A.H. Thompson, Quantitative prediction of permeability in porous rock, Physical Rev. B 34 (1986) 8179-8181.

[40] P.H. Nelson, Permeability-porosity relationships in sedimentary rocks, The Log Analyst 35 (1994) 38-62.

TE D

577

SC

574

[41] S. Ehrenberg, P. Nadeau, Sandstone vs. carbonate petroleum reservoirs: A global

583

perspective on porosity-depth and porosity-permeability relationships, AAPG Bull. 89

584

(2005) 435-45.

586

[42] S.N. Davis, Porosity and permeability of natural materials. In: De Wiest RJM, editor. Flow through porous media, Academic Press, New York (1969) 53-89.

AC C

585

EP

582

587

[43] T. Xu, Y. Ontoy, P. Molling, N. Spycher, M. Parini, K. Pruess, Reactive transport

588

modeling of injection well scaling and acidizing at Tiwi field, Philippines, Geothermics 33

589

(2004) 477-91.

590 591

[44] P. Ungemach, Reinjection of cooled geothermal brines into sandstone reservoirs, Geothermics 32 (2003) 743-761.

592

21

ACCEPTED MANUSCRIPT TABLES Geothermal assessment of the Pisa plain, Italy: coupled thermal and hydraulic modeling Guanhong FENG, Tianfu XU, Fabrizio GHERARDI, Zhenjiao JIANG, Stefano BELLANI

Carbonate reservoir 0.08

Permeability

5.0×10-14 m2

Dry density of the matrix

2700 kg/m3

Specific heat of the matrix

830 J/kg/°C

M AN U

SC

Porosity

Thermal conductivity in reservoir Wellbores Roughness

RI PT

Table 1. Thermal and hydrological parameters of wellbores and carbonate reservoir

3.0 W/m/°C

0.046 mm

3.0 W/m/°C

AC C

EP

TE D

Heat conductivity

1

ACCEPTED MANUSCRIPT FIGURES Geothermal assessment of the Pisa plain, Italy: coupled thermal and hydraulic modeling

EP

TE D

M AN U

SC

RI PT

Guanhong FENG, Tianfu XU, Fabrizio GHERARDI, Zhenjiao JIANG, Stefano BELLANI

AC C

Figure 1. Location map of the area under study in the Pisa plain, with indication of outcropping geological formations and their permeability. The blue contour is the approximate extension of the modelled area.

1

ACCEPTED MANUSCRIPT FIGURES Geothermal assessment of the Pisa plain, Italy: coupled thermal and hydraulic modeling

M AN U

SC

RI PT

Guanhong FENG, Tianfu XU, Fabrizio GHERARDI, Zhenjiao JIANG, Stefano BELLANI

AC C

EP

TE D

Figure 2. Setup of the thermal system and geological framework.

1

ACCEPTED MANUSCRIPT FIGURES Geothermal assessment of the Pisa plain, Italy: coupled thermal and hydraulic modeling

TE D

M AN U

SC

RI PT

Guanhong FENG, Tianfu XU, Fabrizio GHERARDI, Zhenjiao JIANG, Stefano BELLANI

AC C

EP

Figure 3. Lateral 2D cross-section of the domain (a) and 3D view of the model grid (b), together with wellbore structure of the geothermal doublet (c).

1

ACCEPTED MANUSCRIPT FIGURES Geothermal assessment of the Pisa plain, Italy: coupled thermal and hydraulic modeling

TE D

M AN U

SC

RI PT

Guanhong FENG, Tianfu XU, Fabrizio GHERARDI, Zhenjiao JIANG, Stefano BELLANI

AC C

EP

Figure 4. Temperature (a) and pressure (b) distribution after 50 years along a 2D WSW-ENE trending section of the modelled domain (reference case simulation, RCS, with 80 m3/h injection/production constant rate). Red lines are the traces of the wells.

1

ACCEPTED MANUSCRIPT FIGURES Geothermal assessment of the Pisa plain, Italy: coupled thermal and hydraulic modeling

M AN U

SC

RI PT

Guanhong FENG, Tianfu XU, Fabrizio GHERARDI, Zhenjiao JIANG, Stefano BELLANI

AC C

EP

TE D

Figure 5. Temperature distribution after 50 years along a 2D WSW-ENE trending section of the modelled domain (600 m3/hour case). Red dashed lines are the traces of the wells.

1

ACCEPTED MANUSCRIPT FIGURES Geothermal assessment of the Pisa plain, Italy: coupled thermal and hydraulic modeling

M AN U

SC

RI PT

Guanhong FENG, Tianfu XU, Fabrizio GHERARDI, Zhenjiao JIANG, Stefano BELLANI

AC C

EP

TE D

Figure 6. Temperature (a) and heat extraction rate (b) variation for different constant production rates.

1

ACCEPTED MANUSCRIPT FIGURES Geothermal assessment of the Pisa plain, Italy: coupled thermal and hydraulic modeling

M AN U

SC

RI PT

Guanhong FENG, Tianfu XU, Fabrizio GHERARDI, Zhenjiao JIANG, Stefano BELLANI

AC C

EP

TE D

Figure 7. Pressure field after 0.1 and 50 yrs (reference case simulation, RCS, with 80 m3/h injection/production constant rate).

1

ACCEPTED MANUSCRIPT FIGURES Geothermal assessment of the Pisa plain, Italy: coupled thermal and hydraulic modeling

M AN U

SC

RI PT

Guanhong FENG, Tianfu XU, Fabrizio GHERARDI, Zhenjiao JIANG, Stefano BELLANI

AC C

EP

TE D

Figure 8. Time variation of the outflow temperature for two different scenarios. “Scenario 1”: 80 m3/hour constant rate over 12 months per year; “Scenario 2”: 200 m3/hour constant rate over 5 months per year. “Scenario 2” curves plot only data for the 5 month period of production. During the remaining 7 months without production (not shown), temperature goes down to 20°C, and quickly rises to more than 50°C as soon as production re-starts.

1

ACCEPTED MANUSCRIPT FIGURES Geothermal assessment of the Pisa plain, Italy: coupled thermal and hydraulic modeling

M AN U

SC

RI PT

Guanhong FENG, Tianfu XU, Fabrizio GHERARDI, Zhenjiao JIANG, Stefano BELLANI

AC C

EP

TE D

Figure 9. (a) Heat extraction rates predicted by the coupled and uncoupled models in the RCS . (b) Temperature profile along injection and production wells after 50 years (coupled model).

1

ACCEPTED MANUSCRIPT FIGURES Geothermal assessment of the Pisa plain, Italy: coupled thermal and hydraulic modeling

M AN U

SC

RI PT

Guanhong FENG, Tianfu XU, Fabrizio GHERARDI, Zhenjiao JIANG, Stefano BELLANI

AC C

EP

TE D

Figure 10. (a) Lognormal distribution of permeability (variance: 0.8; mean: 5.0 × 10-14 m2. (b) Experimental semivariogram (spherical model) inferred from 100 generated permeability values.

1

ACCEPTED MANUSCRIPT FIGURES Geothermal assessment of the Pisa plain, Italy: coupled thermal and hydraulic modeling

AC C

EP

TE D

M AN U

SC

RI PT

Guanhong FENG, Tianfu XU, Fabrizio GHERARDI, Zhenjiao JIANG, Stefano BELLANI

Figure 11. Three hypothetical heterogeneous permeability distributions (a: HC1, c: HC2, e: HC3) together with the corresponding 2D cross-sections drawn at Y=2500 (b, d, f). Red dashed lines provide the approximate location of the injection and production wells.

1

ACCEPTED MANUSCRIPT FIGURES Geothermal assessment of the Pisa plain, Italy: coupled thermal and hydraulic modeling

TE D

M AN U

SC

RI PT

Guanhong FENG, Tianfu XU, Fabrizio GHERARDI, Zhenjiao JIANG, Stefano BELLANI

AC C

EP

Figure 12. Temperature profile after 50 years in the production well in the reference case simulation (RCS) and in three selected heterogeneous cases (HC1, HC2, HC3).

1

ACCEPTED MANUSCRIPT FIGURES Geothermal assessment of the Pisa plain, Italy: coupled thermal and hydraulic modeling

TE D

M AN U

SC

RI PT

Guanhong FENG, Tianfu XU, Fabrizio GHERARDI, Zhenjiao JIANG, Stefano BELLANI

AC C

EP

Figure 13. 2D WSW-ENE section showing the temperature distribution after 50 years in a homogenous reservoir (a: RCS) and in three hypothetical heterogeneous configurations of the same reservoir (b: HC1; c: HC2; d: HC3). Red dashed lines indicate the approximate location of the injection and production wells.

1

ACCEPTED MANUSCRIPT FIGURES Geothermal assessment of the Pisa plain, Italy: coupled thermal and hydraulic modeling

M AN U

SC

RI PT

Guanhong FENG, Tianfu XU, Fabrizio GHERARDI, Zhenjiao JIANG, Stefano BELLANI

AC C

EP

TE D

Figure 14. Heat extraction (a: variation over the 0 to 50 years interval) and flow (b: snapshot after 50 years) rates in a homogenous reservoir (RCS) and in three hypothetical heterogeneous configurations of the same reservoir (HC1, HC2, HC3).

1

ACCEPTED MANUSCRIPT FIGURES Geothermal assessment of the Pisa plain, Italy: coupled thermal and hydraulic modeling

M AN U

SC

RI PT

Guanhong FENG, Tianfu XU, Fabrizio GHERARDI, Zhenjiao JIANG, Stefano BELLANI

AC C

EP

TE D

Figure 15. Time variation of wellhead pressure in the production (a) and injection (b) well. Four scenarios are compared: homogenous reservoir (RCS) vs. three hypothetical heterogeneous configurations of the same reservoir (HC1; HC2; HC3).

1

ACCEPTED MANUSCRIPT

HIGHLIGHTS “Geothermal assessment of the Pisa plain, Italy: coupled thermal and hydraulic modeling” by Guanhong FENG, Tianfu XU, Fabrizio GHERARDI, Zhenjiao JIANG, Stefano BELLANI.

AC C

EP

TE D

M AN U

SC

RI PT

[1] Feasibility of a low-enthalpy geothermal project assessed by numerical modeling [2] Geothermal district heating project as a part of a local smart city action plan [3] Advanced numerical code applied [4] Coupled groundwater flow and thermal transport [5] Full coupling between wellbore and reservoir