Flow Paths in the Húsmúli Reinjection Zone, Iceland

Examensarbete vid Institutionen för geovetenskaper

Degree Project at the Department of Earth Sciences

ISSN 1650-6553 Nr 428

Flow Paths in the Húsmúli Reinjection Zone, Iceland

Flödesvägar i Húsmúli- återinjektionszonen, Island

Sigrún Tómasdóttir

INSTITUTIONEN FÖR

Examensarbete vid Institutionen för geovetenskaper

Degree Project at the Department of Earth Sciences

ISSN 1650-6553 Nr 428

Flow Paths in the Húsmúli Reinjection Zone, Iceland

Flödesvägar i Húsmúli- återinjektionszonen, Island

Sigrún Tómasdóttir

This thesis work was carried out in cooperation with Reykjavík Energy.

Abstract

Flow Paths in the Húsmúli Reinjection Zone, Iceland Sigrún Tómasdóttir

Reinjection of spent geothermal fluids has become common practise in geothermal power plants. Reinjec- tion can, despite being mostly beneficial, have unwanted effects such as cooling of nearby production wells and injection-induced earthquakes. Tracer tests, along with their modeling and interpretation, are important tools for monitoring the flow paths of the injected water and to predict reservoir cooling.

Knowledge of flow paths in the system allows for better resource management and a more sustainable utilization.

A simulation model of the Húsmúli reinjection zone in the Hellisheiði Geothermal Power Plant in SW-Iceland was developed using the TOUGH2 program. Its hydrological parameters, porosity and permeability, were calibrated using results from an extensive tracer test carried out in the area in 2013- 2015. The aim of the simulations was to obtain better understanding of the flow paths in Húsmúli since, despite fast tracer recovery in production wells in the area, hardly any cooling has been observed in those production wells. The results show that the tracer recovery can be modelled by means of permeable flow channels within the medium. Good results for tracer arrival and concentration peaks were obtained both by assuming a single wide channel and several narrower ones. The parameters that gave the best fit for the single channel model were permeability of 5 · 10 ⁻¹² m ² and porosity ranging from 0.2 % – 3 %. For the multi-channel model they were 1 · 10 ⁻¹² m ² and 0.2 % – 3.5 %, respectively. The high permeability and low porosity in the channels make for an abstract representation of fractured zones within the medium.

Greater cooling was seen with the single-channel modelling approach than with the multiple narrower channel approach, the latter showing hardly any cooling in the production elements during the simulation time. This indicates that the flow paths are more likely multiple channels consisting of fracture networks.

The simulations show that the flow paths are lengthened by sinking of the fluid to greater depth because of the higher density of the colder injected water. This implies that the injected fluid is warmed up by contact with a larger volume of rock, causing a limited and delayed cooling effect.

Keywords: Húsmúli, Hellisheiði Power Plant, geothermal reinjection, tracer tests, numerical modeling Degree Project E1 in Earth Science, 1GV025, 30 credits

Supervisors: Auli Niemi, Gunnar Gunnarsson and Edda Sif Pind Aradóttir

Department of Earth Sciences, Uppsala University, Villavägen 16, SE-752 36 Uppsala (www.geo.uu.se) ISSN 1650-6553, Examensarbete vid Institutionen för geovetenskaper, No. 428, 2018

The whole document is available at www.diva-portal.org

Populärvetenskaplig sammanfattning

Flödesvägar i Húsmúli-återinjektionszonen, Island Sigrún Tómasdóttir

Geotermisk energi anses vara en förnybar och miljövänlig energikälla. Som sådan, kan den spela en viktig roll för att minska utsläppen av växthusgaser från energisektorn över hela världen och genom det bekämpa antropogena klimatförändringar. Geotermiska kraftverk extraherar het vätska från berggrunden, separerar ångan från vätskan och använder sedan ångan för att driva turbiner som genererar elektricitet. Injektion av använd geotermisk vätska från kraftverk har blivit vanligt i den geotermiska industrin för att kassera använd geotermisk vätska, upprätthålla systemtrycket och öka produktionseffektiviteten. Återinjektion av nedkyld vätska kan, trots att den är mestadels fördelaktig, ha oönskade effekter, såsom kylning av närliggande produktionsbrunnar och injektionsinducerad seismisk aktivitet. Spårprov, som möjliggör spårning av en kemikalie inom systemet, tillsammans med modellering, är viktiga verktyg för att förstå flödesvägarna för det injicerade vattnet samt att kunna förutsäga nedkylningar av vattenmagasin. Kunskap om flödesvägar i systemet möjliggör bättre resurshantering och ett mer hållbart användande.

En simuleringsmodell av återinjektionszonen för det geotermiska kraftverket Hellisheiði på sydväs- tra Island, Húsmúli, utvecklades med hjälp av simuleringsprogrammet TOUGH2. Dess hydrologiska parametrar, permeabilitet och porositet, kalibrerades med hjälp av resultat från ett omfattande spårtest som utfördes i området 2013-2015. Syftet med simuleringarna var att få en bättre förståelse av flödesvägarna i Húsmúli. Detta er inressant eftersom trots en snabb återhämtning av spårämne i produktionsbrunnar, har knappt någon kylning observerats i området. Resultaten visar att återhämtningen av spårämnet inte kan modelleras med ett homogent medium, men kan istället modelleras genom att bygga permeabla strömningskanaler inom mediet. Goda resultat för spårämnesankomst och koncentrationstoppar erhölls både genom att använda en enda bred kanal och flera smalare.

Kanalerna ger en abstrakt representation av sprickzoner inom mediet. Större kylning observerades för modelleringsmetoden med en enkel bred kanal än med flera smalare kanaler. Detta indikerar att flödesvägarna i området troligtvis går genom flera sprickzoner. Flödesvägarna förlängs genom att vätskan sjunker till ett större djup på grund av den högre densiteten hos det injicerade vattnet. Detta innebär att den injicerade vätskan värms upp genom kontakt med en större volym berg, vilket medför en begränsad och fördröjd kylningseffekt.

Nyckelord: Húsmúli, Hellisheiði kraftverk, geotermisk injektion, spårtest, numerisk modellering Examensarbete E1 i geovetenskap, 1GV025, 30 hp

Handledare: Auli Niemi, Gunnar Gunnarsson och Edda Sif Pind Aradóttir

Institutionen för geovetenskaper, Uppsala universitet, Villavägen 16, 752 36 Uppsala (www.geo.uu.se)

ISSN 1650-6553, Examensarbete vid Institutionen för geovetenskaper, Nr 428, 2018

List of Figures

1 Schematic geothermal system . . . . 1

2 A map of the Hengill area . . . . 2

3 Rift segments and volcanic zones in Iceland . . . . 6

4 Results of resistivity measurements . . . . 8

5 The resistive core in the Hengill area . . . . 8

6 A map of the Húsmúli reinjection zone . . . . 11

7 Conceptual model . . . . 16

8 Temperature distribution in the system at a depth of 1000 m b.s.l. . . . 16

9 Pressure distribution in the system at a depth of 1000 m b.s.l. . . . 17

10 Illustration of the Voronoi grid generation . . . . 18

11 A map of the numerical grid . . . . 19

12 Rock-type distribution within the layers containing the geothermal system in the model setup . . . . 20

13 Total production rates, injection rates and temperature of injected water during the simu- lation time . . . . 21

14 Location of the feedzones that represent sources/sinks in the model . . . . 23

15 The measured recovery of the 1,3,6-NTS tracer . . . . 25

16 Tracer distribution with the homogeneous system approach . . . . 27

17 Rock-type distribution in the single permeable channel approach . . . . 28

18 Tracer distribution with the single permeable channel approach . . . . 29

19 Tracer recovery with the single permeable channel approach . . . . 29

20 Tracer transport at different times with the single permeable channel approach . . . . 30

21 Best fit for tracer recovery obtained with the single permeable channel approach . . . . . 31

22 Rock-type distribution in the multi-channel approach . . . . 31

23 Tracer distribution with the multi-channel approach . . . . 32

24 Best fit for tracer recovery obtained with the multi-channel approach . . . . 33

25 Temperature distribution at different times in the model . . . . 34

26 Temperature evolution with time in the feedzones in wells HE-31 and HE-48 . . . . 35

27 Temperature profiles at different times for wells HE-31 and HE-48 . . . . 35

Table of Contents

1 Introduction 1

1.1 Geothermal energy and its utilization in Iceland . . . . 1

1.2 Geothermal reinjection . . . . 3

1.3 Reinjection at the Hellisheiði Power Plant, SW-Iceland . . . . 3

2 Background 5 2.1 Site description . . . . 5

2.1.1 Iceland’s tectonic and volcanic setting . . . . 5

2.1.2 The Hengill volcanic system . . . . 7

2.1.3 Study Area: Húsmúli . . . . 9

2.2 Previous research . . . . 11

2.2.1 Numerical modeling in the geothermal industry . . . . 11

2.2.2 Modeling of the Hengill area . . . . 12

2.2.3 Tracer tests in geothermal systems . . . . 13

3 Methods 14 3.1 Theory and software . . . . 14

3.2 Conceptual model development . . . . 15

3.3 Numerical grid generation . . . . 17

3.4 Model setup, boundary and initial conditions . . . . 19

3.5 Simulation Procedure . . . . 23

3.6 Tracer calibration . . . . 24

4 Results 27 4.1 Model calibration . . . . 27

4.1.1 Homogeneous system approach . . . . 27

4.1.2 Single permeable channel approach . . . . 28

4.1.3 Multi-channel approach . . . . 31

4.2 Temperature evolution . . . . 33

5 Discussion 36 5.1 Model approach comparison . . . . 36

5.2 Model assumptions and limitations . . . . 37

5.3 Future studies . . . . 39

6 Conclusions 40

7 Acknowledgements 41

8 References 42

1 Introduction

The aim of this project is to develop a simplified numerical model of the Húsmúli injection zone, and its interaction with the nearby production wells, in the Hellisheiði geothermal field in the Hengill area in SW-Iceland. Furthermore, the aim is to calibrate the hydrological parameters of the model, porosity and permeability, using tracer recovery results from an extensive tracer test that was conducted at the site in June 2013. This kind of model has not been constructed specifically for the Húsmúli formation before.

The main focus will be on injection well HN-17 in Húsmúli and production wells HE-31, HE-44 and HE-48 in Skarðsmýrarfjall mountain which showed the greatest tracer recovery. The aim is to answer the following questions:

• What are the characteristics of the flow paths of the reinjected water in the Húsmúli formation?

• Simple fracture models predict significant cooling in some wells, why do the monitoring data hardly show any cooling?

The software employed is the TOUGH2 simulation program developed at the Lawrence Berkley National Laboratory. The program is widely used for geothermal applications and is used for the overall model of the Hengill geothermal system. The modeling approach will be an effective continuum method where the properties of the medium are chosen so that they roughly represent the fractured porous medium but not the individual fractures themselves.

1.1 Geothermal energy and its utilization in Iceland

Geothermal systems can form where there is a heat source, abundant water and sufficient permeability. In some areas, particularly tectonic or volcanic areas, high temperatures can be found relatively close to the surface which allows for diverse utilization of the geothermal energy. The heat sources in geothermal systems in volcanic areas are cooling bodies of magma in the crust (fig. 1) (Gunnarsson & Aradóttir 2014).

Figure 1. A cross section through a hypothetical geothermal system showing the needed components; a heat source, abundant water and high permeability (Gunnarsson & Aradóttir 2014, Fig. 1).

Geothermal energy is considered a renewable and clean energy source. As such, it can play an important

role in reducing the green house gas emissions from the energy sector worldwide and by that combat

anthropogenic climate change. Geothermal energy has been used directly for applications such as bathing and cleaning clothes since the beginning of civilization but in recent decades, more and more emphasis has been put on electricity generation (Friðleifsson et al. 2008). Power plants extract the hot fluid from the subsurface, separate the steam from the fluid and then use the steam to drive turbines that generate electricity. The power plant operators need to dispose of the separated water and the steam condensate from the turbines. That can be done by reinjecting the fluid back into the ground (Kristjánsson et al. 2016).

A typical geothermal power plant emits only about 5–6 % of the CO ₂ emissions of a fossil fuel power plant of the same size (Aradóttir et al. 2015).

Iceland has vast geothermal resources and harnessing of geothermal energy has played a large role in improving the quality of life in the country. There is a long history of direct use in the country for purposes such as space heating, industry, bathing and even cooking (Gunnlaugsson & Ívarsson 2010). In recent years and decades, there has been increased interest in electricity production and currently (2016), 27.3 % of the electricity production in the country comes from geothermal resources (Orkustofnun 2017).

Seven geothermal power plants are currently operated in Iceland with a combined installed capacity of

663 MWe (Orkustofnun n.d.). ON Power, a subsidiary of Reykjavík Energy, runs the two largest power

plants. They are both located in the Hengill area; in Hellisheiði and Nesjavellir (fig. 2). Drilling of wells in

the Nesjavellir field started already in the mid-sixties and a thermal power plant was commissioned there

in 1990. Because of favorable conditions, interest awoke for an expansion of the exploitation towards

south of the Hengill area, Hellisheiði. Drilling for the Hellisheiði Power Plant started in early 2000s

(Björnsson et al. 2003). The combined electricity generation capacity of both plants is 423 MWe and their

thermal energy production capacity is 420 MWth. Reinjection has been an integral part of the geothermal

utilization in Hellisheiði since the plant was commissioned in 2006 (Kristjánsson et al. 2016).

1.2 Geothermal reinjection

Reinjection of spent geothermal fluids is an important part in all geothermal utilization projects that aim to be sustainable and environmentally friendly. In the beginning, reinjection was simply a way to dispose of spent fluids but with time the purpose of reinjection has changed. Reinjection acts as additional recharge into the system, minimizes the pressure drop in a reservoir subjected to continued utilization and prevents land subsidence. Since most of the thermal energy in geothermal systems (80–90 %) is stored in the reservoir rocks, not the fluid itself, reinjection also serves to extract more energy from the system and by that increase the operational efficiency (Axelsson 2012; Kristjánsson et al. 2016). Reinjection into deep layers can also ensure that injected water containing toxic chemicals does not percolate into freshwater aquifers used to supply drinking water. This can, for instance, be achieved by injecting into the ground below impermeable layers that prevent mixing of the geothermal water and the fresh groundwater (Icelandic National Planning Agency 2006). It can, therefore, be said that reinjection is practiced both for environmental reasons and to improve efficiency in production.

Reinjection can, however, also have less desirable effects. Examples of those are scaling due to mineral precipitation, clogging of aquifers, induced seismicity and cooling of production wells. The risk of cooling can be especially high in fractured reservoirs where the colder injected water can travel through rapid migration paths, such as fractures and faults, and quickly reach production areas (Axelsson 2012;

Bessason et al. 2012). Because of these dangers, geothermal power plant operators strive to thoroughly understand each system by conducting extensive tests and research before reinjection plans are made.

Tracer tests and geothermal modeling play a key role in understanding flow paths in geothermal systems and are valuable for operational purposes and future predictions (Axelsson 2012).

1.3 Reinjection at the Hellisheiði Power Plant, SW-Iceland

The Hellisheiði geothermal field (fig. 2) is a liquid dominated field but boiling takes place due to pressure drop caused by production in the area. The power plant consists of six 45 MWe high pressure and one 33 MWe low pressure turbine generator units. The geothermal fluid is separated into steam, which is used to produce electricity, and water, which is used to heat up cold groundwater. The current thermal energy capacity of the plant is 133 MWth but it will gradually increase to 400 MWth as demand increases in Reykjavík. The separated water and the condensed steam are then reinjected (Gunnarsson et al. 2015). To decrease the emissions of environmentally harmful gases from the power plant, CO 2 and H 2 S emissions from the power plant are captured, dissolved in water and reinjected alongside the geothermal fluid back into the ground (Aradóttir et al. 2015). The temperature of the geothermal fluid upon reinjection is under normal conditions between 60 and 80 ^◦ C. The temperature can, however, reach 120–173 ^◦ C if the low pressure turbine or the thermal plant are under maintenance. Without reinjection and natural recharge, the pressure in the system would drop by more than 10 bar annually because of high production density in the area (Kristjánsson et al. 2016). Reinjection is mainly practised at two sites, Gráuhnúkar and Húsmúli (fig.

2). The Húsmúli site was meant to replace Gráuhnúkar as the main reinjection site since the latter proved

to have temperatures fitting for production. Operation of the Húsmúli site has, however, been somewhat

troublesome. Firstly, because the injection capacity of the zone has been declining, secondly, because the injectivity has turned out to be highly dependent on the temperature of the injected water, and, thirdly, because the injection has triggered seismicity in the area (Gunnarsson et al. 2015).

The complications during the reinjection have raised many questions about the Húsmúli formation.

Extensive tracer tests using naphthalene sulfonate tracers were carried out in 2013–2015, both at Húsmúli

and Gráuhnúkar. The goal with these tests was to study the connection between the reinjection wells

and the production field and to estimate the risk of thermal breakthrough. Fracture models presented

by Kristjánsson et al. (2016) using tracer recovery data from the above mentioned tracer test predicted

significant cooling in some of the wells. Monitoring data has, however, hardly shown any significant

cooling after six years of large-scale injection. The results from the fracture models indicate that the flow

paths are most likely more complex than previously thought and that simple models that do not take the

surrounding medium into account cannot accurately capture the system characteristics. There is a need to

better understand the system to improve the efficiency and functionality of the reinjection system.

2 Background

The first subsections in this chapter will describe the site of the modeling study; from a wider geographic context to a more specific local context. They will also review information about specific wells and describe the execution of the tracer test. The later subsections in the chapter will then review the steps involved in numerical modeling of geothermal reservoirs, present a short overview of modeling in the Hengill area and, lastly, cover the importance of tracer tests in the geothermal industry.

2.1 Site description

2.1.1 Iceland’s tectonic and volcanic setting

Iceland’s complex volcanic and tectonic features are the result of the island’s location astride a divergent plate boundary, the Mid-Atlantic ridge, and on top of a hotspot or a mantle plume. The presence of the hotspot produces a thick crust and a more complicated deformation zone than is observed elsewhere on the ridge (Einarsson 1991, 2008). On land, the ridge is represented by numerous boundary segments that are either purely divergent or oblique and linked together by transform zones. The most notable rift segments are the Western Volcanic Zone (WVZ) and the Northern Volcanic Zone (NVZ). These two zones are offset by the Mid-Iceland Volcanic Zone (MVZ) or the Mid-Iceland Belt (MIB). The Eastern Volcanic zone (EVZ) stretches from the center of the country to the Vestmanna Islands and is propagating to the south (Harðarson et al. 1997). The NVZ is connected to the offshore Kolbeinsey Ridge (KR) by the Tjörnes Fracture Zone (TFZ). The EVZ is then connected to the WVZ by another transform segment, the South Iceland Seismic Zone (SISZ). Lastly, the WVZ is connected to the offshore Reykjanes Ridge (RR) by the Reykjanes Volcanic Belt (RVB), an oblique rift (Harðarson et al. 2009). In the transform segments, volcanism is negligible and strike-slip faulting is most common (Einarsson 2008). Over time, the spreading zone is expected to move completely from the WVZ to the EVZ in a so called ridge jump.

Ridge jumps have been inferred to have taken place earlier in the evolution of Iceland as a result of the

interaction of the mantle plume with the rifting zones (Harðarson et al. 1997). Thirty active volcanic

systems are identified within the volcanic zones in Iceland. They consist of a fissure swarm, a central

volcano or both (Thordarson & Höskuldsson 2008). The different ridge segment, volcanic zones and

volcanic systems can be seen on figure 3.

Figure 3. Location of the rift segments and volcanic zones in Iceland and the distribution of the different volcanic systems. Abbreviations; RR, Reykjanes Ridge; RVB, Reykjanes Volcanic Belt; SISZ, South Iceland Seismic Zone;

WVZ, Western Volcanic Zone; EVZ, Eastern Volcanic Zone; SVB, Snæfellsnes Volcanic Belt; MIB, Mid-Iceland Belt; NVZ, Northern Volcanic Zone; TFZ, Tjörnes Fracture Zone; KR, Kolbeinsey Ridge. The dotted circle shows the approximated center of the mantle plume (Modified from Thordarson & Höskuldsson 2008, Fig. 1).

Extensive geothermal activity is found in Iceland. The geothermal systems are classified into high-

temperature areas, where the temperature at 1 km depth is above 200 ^◦ C, and low-temperature areas,

where the temperature is below 150 ^◦ C (Böðvarsson 1961). The high-temperature areas are found within

the active volcanic zones and the heat source is thought to be a cooling body of magma, either a magma

chamber or dyke swarms (Arnórsson 1995). The low-temperature areas, on the other hand, are mainly

located outside the volcanic zones (Axelsson et al. 2010). The permeability in the high-temperature

systems is highly dependent on the tectonic conditions and seems to be connected to intrusions as well as

faults and fractures (Harðarson et al. 2009). Seismicity within the volcanic zones mostly occurs during

episodes of rifting and magmatism. The largest earthquakes, however, occur on strike-slip faults within

the transform zones (Einarsson 1991).

2.1.2 The Hengill volcanic system

2.1.2.1 Geologic and tectonic setting

Extensive geological exploration has been carried out within the Hengill volcanic system (e.g. Sæmunds- son 1967, 1995a, 1995b; Árnason et al. 1987; Franzson et al. 2005, 2010). The system is composed of a central volcano and a fissure swarm with a graben structure that extends to the northeast and southwest.

The system is about 60 km long and up to 10 km wide (Sæmundsson 2016). Bedrock in the area is primarily composed of hyaloclastite (tuffs, breccias and pillow lavas) formed during glacial periods but also intrusive rock as well as lavas that formed during interglacial periods and flowed to the lowlands (Franzson et al. 2010). The main magma type is tholeiitic basalt although more silicic rock can be found within the Hengill central volcano. The central volcano itself is a hyaloclastite massif composed of tuyas and tindars. Volcanic activity in the system during Holocene is characterized by effusive eruptions on the flanks of the central volcano and on the fissure swarm. The two last eruptions occurred 1900 and 5800 years ago (Sæmundsson 2016). The Hengill system is located at a triple junction between the WVZ, the RVB and the SISZ (fig. 3). This makes for a complicated tectonic setting since the Reykjanes Peninsula is an oblique rift with predominantly NE–SW trending normal faults and volcanic fissures and the SISZ is a fracture zone consisting mostly of N–S trending strike slip faults (Ágústsson et al. 2015). Iceland’s second largest high-temperature geothermal system is located in the Hengill area (Sæmundsson 2016). On the surface, geothermal activity is common as well as extensive thermal alteration that largely correlates with faults and fissures (Franzson et al. 2010).

2.1.2.2 Geophysical exploration

Many geophysical methods have been used to investigate the Hengill geothermal system, e.g. Bouguer

gravity measurements, aeromagnetic measurements, seismic surveys and resistivity measurements. Of

the above mentioned, resistivity measurements are the most informative for delineating the extent of the

geothermal system at depth (Franzson et al. 2010). Extensive MT (Magneto Telluric) and TEM (Transient

Electro Magnetic) studies have been carried out in the area (e.g. Árnason et al. 2010). Geothermal systems

are characterized by a low resistivity cap that marks the outer margins of the reservoir and a more resistive

core within the resource. The resistivity correlates with alteration mineralogy. The low resistivity cap

consists of relatively conductive minerals in the smectite–zeolite zone that correspond to temperatures

between 100 and 220 ^◦ C. At higher temperatures (> 250 ^◦ C) more resistant alteration minerals such as

chlorite and epidote are dominant. The temperature distribution can thus be interpreted from the resistivity

measurements (Árnason et al. 2000). The Hengill low resistivity anomaly that marks the outer margin

of the reservoir covers about 112 km ² (Harðarson et al. 2009). Figure 4 shows resistivity measurement

results at sea level and 500 m below sea level ( m b.s.l.). Figure 5 shows the depth to the top of the resistive

core of the Hengill geothermal system. The figure depicts the top of the resistive core in the Húsmúli area

at roughly 500 m b.s.l. The thickness of the low resistivity layer in this area is 200–250 m (e.g. Magnússon,

Eysteinsson & Árnason 2000; Harðarson et al. 2008).

Figure 4. Results of resistivity measurements (1D inversion of TEM data) at sea level (left) and 500 m b.s.l. (right).

The red hatched areas show the extent of the resistive core (Modified from Árnason et al. 2010, Fig. 8).

Figure 5. Depth to the top of the resistive core in the Hengill area (Modified from Árnason et al. 2010, Fig. 14).

2.1.2.3 Subsurface exploration

Information obtained from drilling has given more concrete data on the subsurface character of the reservoir. Drilling data, drill cuttings, lithological logs, temperature and pressure data, fluid samples and productivity data have, together with surface exploration, drawn up a picture of the geothermal system (Franzson et al. 2010). Down-hole measurements in wells are the only direct measurements of the current primary reservoir parameters. Temperature and pressure data from down-hole measurements are very valuable and can be interpolated to obtain an overview of the whole reservoir. Interpolation between wells will, however, inevitably give rise to uncertainty (Gunnarsson, Arnaldsson & Oddsdóttir 2011).

Temperatures within the resource range from about 200 ^◦ C to about 320 ^◦ C (Franzson et al. 2010). The geothermal system is separated from the colder groundwater system above by the cap rock which has low permeability (Gunnarsson, Arnaldsson & Oddsdóttir 2011). Hydrothermal alteration goes from unaltered rocks in the colder groundwater system, to zeolites in warmer areas and up to high-temperature alteration minerals such as epidote, wollastonite and actinolite within the geothermal system (Helgadóttir et al. 2010). Comparison between estimated formation temperature based on alteration mineralogy and measured temperature logs is used to tell whether the system is in equilibrium or if cooling or heating has occurred (e.g. Árnason et al. 2000; Helgadóttir et al. 2010).

2.1.3 Study Area: Húsmúli

Húsmúli is an interglacial lava shield located on the western flanks of the Hengill volcanic system (fig. 2), on the edge of the Hellisheiði temperature anomaly (Sæmundsson et al. 2016). The series of Húsmúli are composed of hyaloclastites, compound lavas, breccia and postglacial lavas (Khodayar, Axelsson &

Steingrímsson 2015). The cold groundwater table in Húsmúli is at a depth of 30 m but the groundwater table in the deeper geothermal reservoir is at a depth of 200 m (Gunnarsson 2011). This shows that the shallower groundwater system is separated from the geothermal reservoir (Gunnarsson & Aradóttir 2014).

Natural, undisturbed groundwater flow in Húsmúli would roughly be towards the southwest but this flow direction has been overridden because of pressure changes since production started. A pressure high in the injection zone and a pressure low in the production zone have reversed the flow (Khodayar, Axelsson

& Steingrímsson 2015).

Húsmúli has a complicated fracture system. Most notable on the surface are two NNE trending normal faults that run through the eastern part of the formation. The western one runs through the Mógil gully and the eastern one makes up the western part of the Sleggubeinsdalur valley and is called the Húsmúli fault (fig. 6). The appearance of the Húsmúli formation at depth in wells drilled to the east of these faults show that they have a throw to the east of up to 170–200 m (Steingrímsson et al. 1993; Harðarson et al. 2008).

According to the most recent geological map of the area, however, the throw of the Mógil fault has been turned towards the west (Sæmundsson et al. 2016). Reinjection wells drilled into Húsmúli target these faults (Gunnarsson et al. 2015). The largest tectonic features observed at the surface generally do not have the same direction as the faults activated in the seismic events that have occurred recently in conjunction with reinjection into Húsmúli. This indicates that the stress field has changed over time (Bessason et al.

2012). Recent mapping of the fracture pattern in Húsmúli shows a complex fracture network consisting of

NNE trending extensional faults and fault segments as well as shears connected to the SISZ transform zone oriented in various directions (Khodayar, Axelsson & Steingrímsson 2015). This fracture network highly affects the permeability in the Húsmúli which appears to be very anisotropic along the fault lines (Kristjánsson et al. 2016).

2.1.3.1 Tracer test execution and well information

The extensive tracer test conducted during 2013–2015 involved the injection of different napthalene sulfonic acid tracers into six different injection wells. The wells are located in both reinjection sites used for the Hellisheiði Power Plant, Húsmúli and Gráuhnúkar. The tracer injected into well HN-17 in Húsmúli, 1,3,6-NTS, showed by far the greatest recovery and this modeling experiment will solely use 1,3,6-NTS tracer data for the calibration. 100 kg of tracer were dissolved in 4 m ³ of water and injected into HN-17 followed by additional 10 m ³ of water to rinse the tank completely. The injection took place on the 20th of June 2013 and the process took less than 2 hours. Following this, daily sampling started in the nearby production wells and with time included more wells. Sampling continued until the end of 2016 but the sampling frequency decreased with time. The samples were analysed at the Iceland GeoSurvey laboratory by high performance liquid chromatography (Kristjánsson et al. 2016). Detection limits for napthalene sulfonates have shown to be approximately 200 ppt or 0.2 µg/kg using this techology (Rose, Benoit &

Kilbourn 2001).

The production wells that showed the greatest recovery were wells HE-31, HE-48 and HE-44 located NE of the injection zone indicating a principal flow direction parallel to the NNE trending normal faults (Kristjánsson et al. 2016). These wells are all directionally drilled and cased until they enter the geothermal reservoir. Information on these three wells and well HN-17 can be seen in table 1. The location of the wellheads, the horizontal projection of the wellpaths and the identified feedzones in the wells can be seen on figure 6. The feedzones represent permeable structures which can for example be faults, intrusions or intersections between different geological units. The feedzones are classified as large, medium or small based on their relative size. Only large and medium size feedzones are shown on figure 6.

Table 1. Information on injection well HN-17 and production wells HE-31, HE-48 and HE-44. The location of the wellhead is given in the ISN93 coordinate system. Abbreviations; W: wellhead; MD: measured depth; TVD: total vertical depth; DOC: depth of casing.

Coordinates (x,y) W [m a.s.l.] MD [m] DOC [m] TVD [m]

HN-17 383439.34 395731.24 270 2200 636 1780

HE-31 385036.59 396779.94 570 2703 727 2277

HE-48 385059.35 396868.85 570 2288 837 1850

HE-44 385037.10 396863.68 570 2606 837 2340

Figure 6. A map of the Húsmúli reinjection zone and surroundings showing elevation contours, surface fractures, eruptive fissures, production and reinjection wells as well as well paths projected to the surface. Location of large and medium size feedzones in the Húsmúli reinjection wells and the production wells from Skarðsmýrarfjall are also shown. 1:20000 (Data Source: Reykjavík Energy and National Land Survey of Iceland).

The bedrock in the topmost 700 m in well HN-17 consists of basaltic lava layers interbedded with hyaloclastite layers and thinner sedimentary layers. A thick hyaloclastite formation is found from 415 m b.s.l. down to 1454 m b.s.l. where basaltic lavas become dominant again. High temperature alteration starts at a depth of around 600 m b.s.l. with the alteration mineral epidote. Several feedzones are found within the well. The ones below the well casing are a medium one at 512 m b.s.l., a large one at 818 m b.s.l. and a large one at 1429 m b.s.l. (Harðarson et al. 2011). The injected fluid along with the tracer is assumed to exit the well predominantly through the largest feedzones. Wells HE-31, HE-48 and HE-44 are drilled from the same platform on the Skarðsmýrarfjall mountain (fig. 6). The bedrock in these wells is predominantly composed of hyaloclastite until they reach the basement of the Hengill central volcano where basalt becomes dominant. The base is considered to be at a depth of between 900 and 1300 m b.s.l.

(Helgadóttir et al. 2010; Snæbjörnsdóttir 2011).

2.2 Previous research

2.2.1 Numerical modeling in the geothermal industry

Numerical modeling of geothermal reservoirs is a developed technology applied in geothermal research

worldwide. The earliest studies started in the 1970s but the acceptance of the benefit of numerical modeling

in the geothermal industry came with the 1980 Code Comparison Study which tested different geothermal simulators on a set of problems. The technology has since then developed with increased experience, extensive research and increased computational power (O’Sullivan, Pruess & Lippmann 2001). Numerical modeling of geothermal systems consists of solving mass and energy balance equations in a given volume where hydrothermal fluid circulation takes place. It requires the integration of well-established theories of fluid flow and heat transfer with available geological, geophysical and geochemical data. It is thus a multidisciplinary field that requires input from various backgrounds. The main aims of geothermal modeling are to be able to reproduce the conditions in the system, estimate the production capacity and predict the response of the system to production (Gunnarsson, Arnaldsson & Oddsdóttir 2011; Franco

& Vacarro 2014). Before a numerical model can be constructed, a so-called conceptual model of the geothermal reservoir in question needs to be developed. The conceptual model illustrates the structure of the system and the main physical and chemical processes that occur within it. It can for example include;

surface manifestations, flow boundaries, geological features, permeability variations, isotherms, deep inflow zones, zones of boiling etc. Once the conceptual model is ready, the numerical model can be constructed with the aim of representing the conceptual model using mathematical expressions.

The continuous model volume needs to be broken down into smaller volume elements for which calculations of mass and energy balances can be performed on. This is referred to as grid generation.

The elements are then assigned to a specific formation depending on the layering and structures in the conceptual model. Each formation is assigned hydrogeological parameters such as permeability, porosity, density etc. (Pruess 2002). Representing fractured media in a numerical model is a difficult task.

The simplest approach is to assume a porous medium but other approaches include double porosity or representations of a few dominant faults (O’Sullivan, Pruess & Lippmann 2001). Information about fluid properties and their dependence on the thermodynamic changes in the system is needed. Such information is often incorporated into the numerical simulation program used. Appropriate boundary conditions need to be assigned to the edges of the model and information needs to be gathered on sinks and sources in the system. Furthermore, initial conditions for the model at the start of simulation need to be estimated.

This is often the natural state of the geothermal system before production. To verify the accuracy of the model, it needs to be calibrated using all available data, such as injection and production rates, pressure, temperature and tracer data (Pruess 2002). The calibration is an iterative process where the aim is to match the model results with real data. This is a challenging task requiring multiple model runs and adjustments to the model parameters (Burnell et al. 2012).

2.2.2 Modeling of the Hengill area

In conjunction with drilling, extensive geological, geophysical and geochemical surveys have been done

in the Hengill area, leading to the existence of a comprehensive data collection. Reservoir modeling has

Hellisheiði shed new light on the geothermal activity in the area and challenged the idea of a common heat source. This led to a revision of the conceptual model and subsequently updates to the numerical model.

The current model assumes that various local heat sources drive the system and that the heat sources are cooling magmatic intrusions and a deeper magma chamber. The model is calibrated using geological survey data, down-hole measurements and production history (Gunnarsson, Arnaldsson & Oddsdóttir 2011). All reservoir models for the Hengill area have been developed using the TOUGH2 family of codes (Björnsson et al. 2003).

2.2.3 Tracer tests in geothermal systems

A tracer is a chemical compound that can be injected into a hydrological system and used to better understand the properties of the system. Its recovery with time at different observation points around the injection well is monitored and used to delineate the flow of water in the system (Axelsson 2012).

Tracer testing has become an important part of geothermal field research and management, especially in reinjection studies (e.g. Rose et al. 1997; Shook 2001; Axelsson 2012). This is because tracer tests can give information about the nature of the flow paths between injection and production zones. For example whether there exist rapid migration paths, such as fractures and faults, between the zones. This is important to assess the risk of cooling from the injection of colder fluid into a geothermal system. Tracer transport is significantly faster than the movement of a cooling front, so tracer test results can be used in a predictive manner (Axelsson 2013). The tracers used depend on the properties of each system. In general, tracers used in geothermal systems should be environmentally friendly, non-reactive, heat tolerant, inexpensive and easy to detect. Furthermore, the tracers should not be present in the system beforehand, or at least in very small quantities. Examples of tracers commonly used in geothermal systems are fluorescein, benzoic acid and napthlene sulfonates (Axelsson 2012). Napthalene sulfonates have proven to be excellent tracers for high temperature geothermal systems since they are thermally stable in temperatures exceeding 310 ^◦ C (Rose, Benoit & Kilbourn 2001).

The recovery of the tracer can be used to qualitatively evaluate the connection between the injection and production wells. That kind of information can then be used in decision making for new injection well locations. Quantitative and comprehensive interpretation of tracer tests is a more complicated task.

Axelsson, Björnsson and Montalvo (2005) presented a method for such interpretation that assumes the existence of specific flow channels that connect the injection and production well feedzones. Results from the method can be used to estimate the characteristics of the flow channel. This method was used to analyze the 2013–2015 tracer test in Hellisheiði and the results were presented by Kristjánsson et al.

(2016). The method showed a good correspondence between the measurements and the simulated values.

However, the results showed large dispersion in the channels indicating fracture networks rather than

simple fracture flow paths. The cooling predictions made using these models predicted cooling of between

10 and 20 ^◦ C after 5 years of injection in wells HE-31 and HE-48 (Kristjánsson et al. 2016). Monitoring

data has, as previously mentioned, hardly shown any cooling since large scale injection started in Húsmúli

in late 2011. The latest measurements in well HE-31 show a cooling of 0.7 ^◦ C in this well between 2015

and 2017 (Gunnarsson 2017).

3 Methods

The following subsections describe the methods used to construct the model and run the flow and tracer simulations. The first section introduces the governing equations and the software used. The second section describes the conceptual model development. The following chapters describe the numerical model setup and simulation strategy.

3.1 Theory and software

In this modeling work, the TOUGH2 simulator (Pruess, Oldenburg & Moridis 2012), as implemented in forward mode in the iTOUGH2 program, was used to iteratively solve the governing equations forward in time (Finsterle 2007). TOUGH2 is a multiphase flow and transport simulation program for fractured and porous media. The governing equations solved by TOUGH2 describe the conservation of mass and energy.

The change in mass/energy in a given subdomain V n resulting from fluxes across enclosing surface Γ n is represented as:

d dt

Z

V

M ^κ dV n = Z

Γ

F ^κ • n dΓ n + Z

V

q ^κ dV n (3.1)

where M ^κ stands for the mass/energy of the mass/heat component κ present in that subdomain. F ^κ • n dΓ n

stands for the flux of component κ into domain V n normal to surface Γ n . Lastly, q ^κ stands for sinks or sources of component κ in domain V n (Pruess, Oldenburg & Moridis 2012).

In the program, conduction and convection control the heat flow. Thermodynamic conditions are based on local equilibrium of all phases. Advection controls the mass flow and a multiphase version of Darcy’s law is used to calculate advective mass fluxes in each phase (equation 3.2).

F _β = ρ _β u _β = −k k _rβ ρ _β

µ _β (5P _β − ρ _β g) (3.2)

where F _β is the flux of phase β, ρ _β is the density of phase β, u _β is the Darcy velocity in phase β, k is the absolute permeability, k _rβ is the relative permeability of phase β, µ _β is the viscosity of phase β, P _β is the sum of the pressure of a reference phase and the capillary pressure and g is the vector of gravitational acceleration. For further specification of the equations solved by TOUGH2, the reader is referred to the manual (Pruess, Oldenburg & Moridis 2012).

Standard TOUGH2 does not take dispersion into account and in this work, diffusion was not taken into account either. However, it should be noted that the numerical grid itself creates some dispersion.

Dispersion modeling in TOUGH2 is possible with a specific module but that is limited to two-dimensional

flow problems. Discretization of the continuous space is made using the "integral finite difference" method.

allows for separate tracking of two water components which contain different trace constituents but have identical physical properties. The second water component can be considered a tracer and its transport though the system can be monitored. Vapor pressure lowering from capillary and adsorption effects is neglected in EOS1 meaning that vapor pressure in two-phase conditions is equal to the saturated vapor pressure of the bulk liquid (Pruess, Oldenburg & Moridis 2012). As previously mentioned, the Húsmúli system is liquid dominated. The primary variables required to define the system in EOS1 are pressure and temperature for single-phase conditions and pressure and gas saturation for two-phase conditions. If two water components are included, the mass fraction of the second water component in each element needs to be defined as well.

3.2 Conceptual model development

The conceptual model developed for this project is a simplified version of the Húsmúli injection site (fig.

7). A 200 m thick impermeable cap rock is assumed to separate a colder groundwater system from the geothermal system causing very little interaction between the two. The boundary between the cap rock and the geothermal system is assumed to be at a depth of 500 m b.s.l. everywhere in the model. The depth and thickness are based on resistivity measurements which show the top of the resistive core in Húsmúli at around 500 m b.s.l. (fig. 5) (Árnason et al. 2010). Alteration mineralogy in well HN-17 also shows the incoming of epidote at 612 m b.s.l. (Harðarson et al. 2011) which indicates temperatures above 230 ^◦ C (Helgadottir et al. 2010). Injection of fluid takes place in Húsmúli into wells HN-14, HN-09, HN-12, HN-16 and HN-17 (fig. 6). Information on the location and depth of feedzones in all wells was supplied by Reykjavík Energy (fig. 6). As a simplification, the fluid is assumed to enter the system predominantly through the largest feedzones in each well. The tracer is injected into well HN-17. Fluid is produced from the wells drilled in Skarðsmýrarfjall, HE-31, HE-48, HE-44 and HE-33 (fig. 6), and the production is also assumed to come from the largest feedzones in those wells.

Temperature and pressure distribution in the system at the beginning of the simulation was extracted

from the existent large scale Hengill model and provided by Reykjavík Energy. This model is based on

all available data from the geothermal field. Figures 8 and 9 show the distribution of temperature and

pressure at a depth of 1000 m b.s.l. for the 15th of September 2011. This point in time marks the start of

large-scale injection in Húsmúli and will be the starting point for the simulations in this study.

Figure 7. A cross-section in the direction of flow showing injection well HN-17 and the three wells that showed

the best recovery, HE-31, HE-48 and HE-44. The wellpaths have been projected to a vertical plane. All identified

medium and large size feedzones in these wells are shown as well as the layer division in the numerical model. The

location of the cross-section is shown on figure 6.

Figure 9. Pressure distribution in the system at a depth of 1000 m b.s.l. The pressure values are interpolated between the point data obtained from the overall Hengill model using the software Paraview. 1:25000 (Data Source:

Reykjavík Energy).

The Húsmúli formation is on the edge of the high-temperature anomaly and no specific influx of heat is assumed in this area of the model during the simulation time period, only a constant bottom temperature.

Húsmúli lies on the edge of the fracture zone that extends through the Hengill area (fig. 2). The part of the model that lies within the fracture zone is assumed to have higher permeability than the surrounding formation. Well HN-11 was meant to be a reinjection well but proved to have very low permeability (Gunnarsson et al. 2016). For this reason, well HN-11 is considered outside the high permeability zone.

Tracer recovery from the 2013–2015 tracer test was very limited in production well HE-46 and wells to the southeast of it (Kristjánsson et al. 2016). Because of this and contrasts in the temperature and pressure distribution, a permeability barrier is assumed between the Húsmúli area and the production zone to the southeast of it.

3.3 Numerical grid generation

The program AMESH was used for the numerical grid creation (Haukwa 1998). It generates numerical

grids for modeling of flow and transport in programs with formulation based on the integral finite difference

method. The Voronoi or Thiessen tessellation method is used to create the elements. The principle behind

the method is that the interfaces between adjacent elements are the perpendicular bisectors of the line

between two points. Element vertices form where the different bisectors intersect. The shape of each

element is thus defined by the points that are closest to each element’s center point (fig. 10).

(a) (b) (c)

Figure 10. Illustration of the Voronoi grid generation (a) Points (b) Lines connecting the points (c) Perpendicular bisectors of lines form elements.

The input into the program is a list of center points, a defined outer boundary of the grid and a connection

length tolerance setting. Each center point has defined x,y,z coordinates as well as element thickness. 3D

grids have vertical interface areas that are horizontal projections of the 2D grid (Haukwa 1998). That

is, the stack of layers in the grid is cut vertically with the shape of each element forming elements with

identical top and bottom areas but possibly varying height depending on the layer thicknesses. The center

points were constructed manually using ArcGIS. The grid was refined at the locations of the feedzones

in well HN-17 to allow for accurate calculation of tracer concentration since it will be highest closest to

the injection points. Mesh refinement around areas of mass transfer is commonly practised in numerical

modeling (Franco & Vacarro 2014). The grid was made coarser away from the injection well in the

direction of flow towards the production zone. The elements at the outer edges of the grid were given

volumes that are orders of magnitude larger than the elements within the active flow zone. Each layer in

the model consists of 1738 elements. The model has 22 layers, labeled from A to V (fig. 7), and thus

consists of 38236 elements. AMESH assigns each element a name according to the TOUGH2 naming

convention. The names are 5 digits and the first one represents the layer the element belongs to. The

layers in the model have varying thickness. Layers that contain feedzones in well HN-17 were made

thinner than others and layers considered of less importance in this study, such as the colder groundwater

system and the bottom layer, were made thicker. Figure 11 shows the distribution of the elements within

the center of the grid and its location in SW-Iceland.

Figure 11. Distribution of the elements within the center of the grid and its location in SW-Iceland. 1:25000 (Data source: Reykjavík Energy and National Land Survey of Iceland).

3.4 Model setup, boundary and initial conditions

Limited data exists below 2500 m b.s.l. and in addition to that the temperature and pressure conditions at

this depth cause the fluid to be in supercritical condition which falls out of the scope of the numerical

simulator used (Gunnarsson, Arnaldsson & Oddsdóttir 2011). For this reason the model only reaches

down to 2500 m below sea level. The model stratification can be seen on the cross-section in figure

7. The geothermal system extends through layers D to U. The model elements were assigned different

rock types using the TOUGH2 preprocessor Steinar (Vatnaskil, 2015). All rock types were given a

thermal conductivity of 2.1 W/(mK) (Clauser & Huenges 1995), heat capacity of 1·10 ³ J/(kgK) (Bouhifd

et al. 2007), density of 2650 kg/m ³ and porosity of 10 % (Gunnarsson & Aradóttir 2014). The initial

permeability values for the general model setup are shown in table 2. The initial permeability in the

geothermal system is chosen as the permeability of hyaloclastite as presented by Aradóttir et al. (2012)

since the system is predominantly composed of hyaloclastite as described in subsection 2.1.3.1.

Table 2. Permeability (k) for the various rock formations defined in the model. The values are obtained from Gunnarsson and Aradóttir (2014) and Aradóttir et al. (2012) except for the feedzones and the permeability barrier which are the author’s initial guess.

Rock Type k [m ² ]

Bottom Layer 1 x 10 ⁻¹⁷

Surrounding formation 5 x 10 ⁻¹⁷ Geothermal system 3 x 10 ⁻¹⁴

Cap rock 1 x 10 ⁻¹⁷

Colder groundwater system 1 x 10 ⁻¹⁴

Top layer 1 x 10 ⁻¹⁶

Feedzones 1 x 10 ⁻¹²

Permeability Barrier 5 x 10 ⁻¹⁵

The top layer, the colder groundwater layer, the cap rock layer and the bottom layer are assigned the values listed in table 2 over their entire horizontal extent. The layers covering the geothermal system as depicted in figure 7 are given a higher permeability within the the active fracture zone but a lower permeability in the surrounding formation and in the permeability barrier. Figure 12 shows an image obtained from the Steinar program of the initial rock-type distribution within the geothermal system. The part of the geothermal system that lies within the Húsmúli area is given another rock type than the part that lies on the other side of the permeability barrier to allow for modifications of parameters solely within Húsmúli.

Figure 12. Rock-type distribution within the layers containing the geothermal system in the model setup. Light

conditions. TOUGH2 offers the possibility of making such elements or layers behave as inactive during a simulation. No energy or mass balance equations are set up for inactive elements and their thermodynamic conditions remain unchanged during the simulation period (Pruess, Oldenburg & Moridis 2012).

The temperature and pressure values at the start of the simulation for the whole model were obtained from the overall model for the Hengill area as previously mentioned. These values serve as initial conditions for the model. The overall Hengill model covers a larger area and has a coarser numerical grid than the model used in this study. Because of this, the data from the Hengill model needed to be manipulated to get values for the elements in this study. Firstly, the data for the elements that lie within the model boundary of this study were extracted. Secondly, an optimal Delaunay triangulation was performed and the data was gridded to a 25x25 m grid. Thirdly, values at the element center coordinates in this model were extracted from the grid. Lastly, the data were linearly interpolated with regards to depth to be able to extract values for the depths of the layers in this model. This procedure was performed for values taken out of the overall Hengill model as it simulates the system conditions on the 15th of September 2011.

Injection history from all active reinjecton wells in the Húsmúli area is incorporated in the model, since it is the total amount of injected and produced fluid that in the end drives the flow in the system. The best recovery of the 1,3,6-NTS tracer injected into HN-17 was in production wells HE-31, HE-48, HE-44 and HE-33 (fig. 6). Production history from these wells is incorporated in the model. Monthly mean values of injection and production rates were supplied by Reykjavík Energy as well as the temperature (or enthalpy) of the injected fluid. The injection and production rates over the simulation time as well as the temperature of the injected water is shown on figure 13.

Figure 13. Total production rates, injection rates and temperature of injected water during the simulation time.

The wells as such are not defined in TOUGH2. Instead, injection or production is defined as sources/sinks in the model elements that coincide with the location of the largest feedzones in each well. The simplifi- cation made in this modeling study is to divide the total flow rate evenly between the largest feedzones in each well. Exceptions to this rule were made for wells HE-44 and HN-12 which only have one large feedzone defined and well HE-33 which only has medium size feedzones defined. In wells HE-44 and HN-12, one medium feedzone was incorporated as well and assigned 20 % of the production/injection rate. In well HE-33, the total production rate was divided between the medium size feedzones within the geothermal system. Elements that represent sources or sinks are assigned a specific feedzone rock type with higher permeability since they represent permeable structures. To account for varying element size in the model, especially between the production elements, which have a radius of about 50 m, and the feedzone elements in well HN-17, the feedzone rock type was assigned to elements in a 50 m radius from the injection elements. Information about the sources and sinks defined in the model is presented in table 3. The sources/sinks are annotated with the well name preceded by a letter running through the sequence A,B,C,D depending on the number of feedzones, from the shallowest to the deepest. The location of the feedzones that represent sources/sinks in the model is shown on figure 14.

Table 3. Annotation, element name, depth and total flow ratio defined for the sources and sinks in the model.

The sources/sinks are annotated with the well name preceded by a letter running through the sequence A,B,C,D depending on the number of feedzones, from the shallowest to the deepest. The element names are assigned by AMESH during the grid generation according to the TOUGH2 naming convention.

Annotation Element name Depth [m b.s.l.] Flow ratio

Sources AHN-09 TB928 1908 1/2

BHN-09 UB572 2091 1/2

AHN-12 EB930 545 1/5

BHN-12 QA883 1466 4/5

AHN-14 KB956 997 1/2

BHN-14 PC167 1419 1/2

AHN-16 OB537 1366-1384 1/4

BHN-16 QB524 1480 1/4

CHN-16 RB524 1524 1/4

DHN-16 RB491 1603 1/4

AHN-17 IA686 818 1/2

BHN-17 PB189 1429 1/2

Sinks AHE-31 LB735 1007-1047 1/3

BHE-31 MB722 1205 1/3

CHE-31 RB819 1582 1/3

AHE-33 EB849 552-572 1/2

BHE-33 HB859 741 1/2

AHE-44 EB810 545 1/5

BHE-44 OB821 1401 4/5

Figure 14. Illustration of the feedzones (blue squares) that represent sources/sinks in the model. The wells and feedzones are shown within the center of the grid. Injection wells are shown in blue and production wells in red.

The initial pressure distribution in the model is shown with transparency.

3.5 Simulation Procedure

The simulation is run in three stages. The first stage is a single water component flow simulation from the 15th of September 2011 until midnight on the 20th of June 2013. This first stage uses the initial conditions described above as well as the injection and production data for all wells shown on figure 14. The aim of this simulation is to obtain initial conditions at the time of the tracer injection. The tracer injection is made by introducing a second water component which consists solely of tracer. The output from the first stage simulation is used to define the primary variables for each element at the start of the second stage simulation. These primary variables are temperature and pressure, but in the case of a two-component simulation, the mass fraction of the second water component needs to be defined for each element as well.

To avoid convergence problems resulting from zero values, a small initial mass fraction of 1 · 10 ⁻¹⁰ kg/kg is defined for the second water component in each element except for the inactive edges.

The second stage is a two-component simulation that simulates the tracer injection. This injection

of the second water component lasts for two hours, which is the time the tracer injection process took

in reality. In reality, 14 m ³ of water were injected in addition to the normal geothermal fluid during

these two hours (approx. 2 kg/s). In the model, this additional amount of water is considered negligible

compared to the 92.73 kg/s of geothermal fluid injection that was taking place in HN-17 during this

period. In the second stage simulation, the normal geothermal fluid injection into the two large feedzones

in well HN-17 (fig. 14) is simply replaced by an injection of solely tracer (water component 2). This

results in a larger amount being injected than the 100 kg that were injected in reality but the results will be

scaled to account for the difference (see chapter 3.6). Injecting a larger amount of tracer facilitates the numerical calculations since small concentrations can cause instability in the simulations. Geothermal fluid injection and production takes place in the other wells in the system during this simulation. To maintain a steady tracer background value and prevent heavy dilution in all injection elements, which causes severe convergence problems, a very small amount of tracer (1 · 10 ⁻¹⁰ times the normal injection) is coinjected with the geothermal fluid in the injection elements in wells HN-09, HN-12, HN-14 and HN-16 during this simulation.

The output from the tracer injection simulation is used as the initial condition for the last stage. The last stage is a two-component simulation from the 20th of June until the 15th of September 2015 where the flow in the system transports the tracer around. In this simulation the injection into HN-17 is changed back into a normal geothermal fluid injection. The small coinjection of tracer mentioned above is continued in this simulation, now in well HN-17 as well. Outputs from the third stage simulation at different times are used to monitor the transport and recovery of the tracer in the different production wells. Data for plotting is extracted from the output from TOUGH2 using the EXT program available at the TOUGH official website. The data are then plotted using the Paraview open-source visualization software.

3.6 Tracer calibration

Figure 15 shows the observed recovery of the 1,3,6-NTS tracer as presented by Kristjánsson et al. (2016).

The tracer arrived in well HE-31 after 14 days and showed the largest concentration peak there. Only 3 days later it arrived in well HE-48 and 35 days after that in well HE-44. The arrival in well HE-33 came 62 days later and was very limited. The fast arrival of the tracer strongly indicates heterogeneity and anisotropy in the system. This recovery data will be used to calibrate the simulation model.

The concentration of the second water component in the TOUGH2 output file is given as a mass fraction (kg/kg) but the tracer recovery data are reported in µg/kg. To account for this and the larger injected amount of tracer, conversions are made before plotting to facilitate comparison of simulation results with observations. In the model, the tracer is injected for 7000 s and during that time, the combined injection rate into well HN-17 is 92.73 kg/s. That gives a total injected tracer amount of:

M = t · Q = 7000 s · 92.73 kg/s = 649110 kg (3.3) As previously explained, only 100 kg of tracer were injected in reality which needs to be taken into account when interpreting the results. That gives a ratio of:

100kg

649110kg = 1.54 · 10 ⁻⁴ (3.4)

Using this ratio, the concentration reported in the TOUGH2 output file is converted into a tracer concen-

The measured recovery shown on figure 15 is from samples from the well head and not from the individual feedzones. To account for this, a weighted average of the concentration in each feedzone with regards to the flow rate defined from that same feedzone is calculated. This is done for each production well.

Equation 3.6 shows an example calculation for a well with three feedzones.

C total = C A · Q _A +C B · Q _B +C C · Q _C Q _A + Q B + Q C

(3.6) where C represents concentration in µg/kg and Q represents flow rate in kg/s.

Figure 15. The measured recovery of the 1,3,6-NTS tracer in the production wells on Skarðsmýrarfjall and to the east of Húsmúli (Kristjánsson et al. 2016, Fig. 3).

Different approaches will be used to calibrate the model. The first approach tests whether the induced

pressure gradient from the injection and production is sufficient to channel the flow of the tracer towards

the production wells and produce a recovery curve. This approach will simply use the initial setup as

shown in figure 12. The second approach tests whether a more permeable flow channel that extends

over the horizontal and vertical extent of the main feedzones in wells HN-17, HE-31, HE-48, HE-44

and HE-33, results in a more realistic recovery. The third approach tests whether making multiple more

permeable flow channels that directly connect the individual feedzones in the wells results in a more

realistic recovery. The permeability and porosity in the channels are manually adjusted with the aim of

attaining a satisfactory match. The permeability is gradually increased from the initial value presented for

the geothermal system in table 2 and the porosity is decreased until the curves for the observed data and

the simulated data show a satisfactory match.

The flow channels are not meant to represent specific fractures but rather a fractured zone that the fluid can quickly travel through. Fault zones are generally composed of a fault core, a damaged zone and the protolith. The fault core is then the part where most displacement takes place and the damaged zone is the network of structures around it which may enhance the fault zone permeability (Caine, Evans and Forster 1996). The idea behind the channel approach in this model is that fault zones that are aligned more or less parallel to the direction of flow in a system can act as flow conduits but also as barriers to flow perpendicular to the fault zone alignment (Caine, Evans & Forster 1996; Khodayar, Axelsson

& Steingrímsson 2015). To represent the lower permeability across the fault zone than along the fault

zone, the permeability at the outer margins of the channels is assigned the same value as the surrounding

formations but not as high as the permeability within the channel. This is defined by assigning a different

permeability value to the element connections at the channel edges in the TOUGH2 input file. The size of

the elements in the model controls the possible width of the flow channels and by that the cross-sectional

area of the fracture zone. To simulate flow through a smaller cross sectional area, as would be the case in

fracture dominated flow, the porosity in the channel can be decreased.