# Geothermal potential of small sub-volcanic intrusions in a typical Icelandic caldera setting

Steffi Burchardt, Mohsen Bazargan, Einar Bessi Gestsson, Christoph Hieronymus, Erika Ronchin, Hugh Tuffen, Michael Heap, Jonathan Davidson, Ben Kennedy, Alex Hobé, et al.

► To cite this version:

Steffi Burchardt, Mohsen Bazargan, Einar Bessi Gestsson, Christoph Hieronymus, Erika Ronchin, et al.. Geothermal potential of small sub-volcanic intrusions in a typical Icelandic caldera setting. Volcanica, 2022, 5 (2), pp.477-507. 10.30909/vol.05.02.477507. hal-03915051

HAL Id: hal-03915051

https://hal.science/hal-03915051

Submitted on 2 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Geothermal potential of small sub-volcanic intrusions in a typical Icelandic caldera setting

Steffi Burchardt* α, Mohsen Bazarganα,β,γ,δ, Einar Bessi Gestssonα,δ, Christoph Hieronymusα,

Erika Ronchinα,ε, Hugh Tuffenζ, Michael J. Heapη,θ, Jonathan Davidson¹, Ben M. Kennedy¹,

Alex Hobéα, and Elodie Saubin¹

α Department of Earth Sciences, Uppsala University, Uppsala, Sweden.

β Geosigma, rock section, Sankt Eriksgatan 113, 113 43 Stockholm, Sweden.

γ Institute for Geotechnical Engineering, ETH Zurich, Zurich 8093, Switzerland.

δ WSL Institute for Snow and Avalanche Research SLF, 7260 Davos Dorf, Switzerland.

ε Icelandic Meteorological Office, Reykjavik, Iceland.

ζ Department of Earth Sciences, Sapienza University of Rome, P.le Aldo Moro 5, 00185 Rome, Italy.

η Lancaster Environment Centre, University of Lancaster, United Kingdom.

θ Université de Strasbourg, CNRS, Institut Terre et Environnement de Strasbourg, UMR 7063, 5 rue René Descartes, Strasbourg F-67084, France.

l Institut Universitaire de France (IUF), Paris, France.

k School of Earth and Environment, University of Canterbury, Christchurch, New Zealand.

## ABSTRACT

Geothermal exploration targets large magmatic intrusions as heat sources because of their size, longevity, and amount of stored energy, but as shallow volcanic plumbing systems comprise numerous smaller intrusions, their geothermal potential warrants consideration. Here, we evaluate the geothermal impact of dykes and sills on caldera-infill rocks. We present geological data and geothermometry on intrusions in the eroded Breiðuvík caldera in Northeast Iceland, which serves as an analogue to the active, and geothermally exploited, Krafla volcano. These data inform 2D finite element models of dyke and sill intrusions that consider heat transfer in porous media. Our results indicate that small intrusions create considerable thermal anomalies in their immediate vicinity. These anomalies are larger-magnitude and longer-lasting for individual thick sills and dykes, but networks of smaller sills and dykes emplaced close in time and space can create more widespread thermal anomalies that may be viable economic targets for decades after their emplacement.

## ÁGRIP

Jarðhitarannsóknir miða oftast að því að finna stór kvikuinnskot sem hitagjafa í jarðhitakerfum vegna lengd storknunartíma og magni af orku. Grynni hluti kvikuflutningskerfa eldfjalla samanstendur af fjölda minni kvikuinnskota sem gefur tilefni til frekari rannsókna á þeim sem hitagjafa. Þessi rannsókn metur áhrif lóðréttra (gangar) og láréttra (sillur) kvikuinnskota á jarðhitakerfi í öskjum. Jarðfræði og efnasamsetning kvikuinnskota í Breiðuvíkur eldstöðinni, sem er rofin megineldstöð staðsett á norðausturhluta Íslands, er kynnt og notuð sem hliðstæða við jarðhitakerfið í Kröflu sem í dag er nýtt til orkuframleiðslu. Gögnin eru notuð í tvívíð reiknilíkön af kvikuinnskotum sem áætla varmaflutning í gropnu bergi. Niðurstöður rannsóknarinnar benda til þess að minni innskot valdi umtalsverðu hitafráviki í sínu næsta umhverfi. Slík frávik eru hærri og langlífari fyrir einstaka þykkari innskot. Kerfi af fjölda minni innskota mynduð nærri í tíma og rúmi geta framkallað umfangsmeiri hitafrávik sem gætu verið raunhæf til jarðhitavinnslu áratugum eftir myndun innskotana.

KEYWORDS: Geothermal energy; Magmatic intrusion; Sill; Dyke; Caldera volcano; Numerical modelling.

## 1 INTRODUCTION

In Iceland, where geothermal power is the main source of primary energy used by society [Orkustofnun - National Energy Authority of Iceland 2021], caldera volcanoes and other volcanic structures are exploited for geothermal energy. Alongside widespread use of geothermal water for heating, electricity is produced in geothermal power plants located in high-temperature geothermal fields, where heat originates from upper crustal magmatic intrusions [Arnórsson et al. 2008]. Currently, exploration and exploitation of supercritical geothermal resources above shallow level magmatic intrusions has the potential to become an even more efficient way of energy pro-

duction [Scott et al. 2015]. Geothermal and scientific drilling in high-temperature geothermal areas testifies to the existence of both large (km-scale) and small (metre- to tens-of-metres scale) magmatic intrusions [Fridleifsson and Elders 2005; Reinsch et al. 2017]. The setting of these intrusions is usually complex. The intrusions are often hosted within caldera volcanoes that have seen multiple stages of activity, including the construction of an eruptible volcanic plumbing system, characterised by eruptions that vary in size and style, and experience (partial) destruction/collapse of the volcanic edifice and the underlying plumbing system during caldera formation, as well as rejuvenation with renewed intrusive and extrusive activity [Kennedy et al. 2018]. As a result, calderas are complex patchworks of different rock types, structures, and intrusions.

<figure><img src="image_2.png" /></figure>

steffi.burchardt@geo.uu.se
Typical caldera-filling rock types may comprise lava flows, pyroclastic deposits, lacustrine sediments, and even subglacial deposits. Needless to say, the physical properties of these diverse rock types, such as porosity and permeability, differ considerably, as revealed in previous laboratory studies [e.g. Farquharson et al. 2015; Schaefer et al. 2015; Heap et al. 2017; 2020c]. Furthermore, these rocks are also affected by caldera collapse-related and/or tectonic faulting [e.g. Burchardt and Walter 2010], and are invaded by magmatic intrusions, driving hydrothermal alteration and mineralisation [e.g. Morden sky et al. 2018; Heap et al. 2020a; Kennedy et al. 2022]. Typical intrusion types in calderas include cryptodomes/accoliths, cauldron-subsidence plutons, and magmatic sheet intrusions, such as dykes and sills [Kennedy et al. 2018]. As a result of the abundant magmatic activity, and the variety of structural and lithological barriers and reservoirs that favour the collection and availability of hot fluids, calderas are ideal sites for mineral and geothermal exploration [e.g. Wohletz and Heiken 1992; Stix et al. 2003; Garden et al. 2017] with high temperatures possible in a variety of stratigraphic and structural settings [Garden et al. 2020].

Geothermal areas in Iceland have been modelled numerically to simulate the behaviour of the systems during geothermal production [e.g. Björnsson 1999; Björnsson et al. 2003; Gunnarsson et al. 2011]. However, such models frequently exclude the heat source of the geothermal field, the solidifying intrusions themselves. These intrusions are instead assumed to be out of the depth range of the models and the heat source compensated by adjusting the boundary conditions of the bottom layer (heating the entire model bottom) [e.g. Karaoğlu et al. 2019]. Although this assumption often works for practical purposes, in reality, magmatic intrusions reach well into the depth range of geothermal fields, as evident from events where geothermal wells were either used as magma pathways (e.g. during the Krafla Fires) [Larsen et al. 1979] or when comparatively shallow boreholes intersected molten intrusions [Elders et al. 2011; 2014]. In studies where the heat source is included in numerical models, fluid flow around magma bodies mostly addresses pluton geometries on the scale of several kilometres [e.g. Hayba and Ingebritsen 1997; Eldursi et al. 2009; Annen 2017], which is not representative for shallow geothermal systems where the majority of intrusions are far smaller-volume dykes and sills [cf. e.g. Kennedy et al. 2018]. Hence, geothermal modelling needs to account for the local nature of shallow magmatic heat sources [Gunnarsson and Aradóttir 2015].

This study is based on a two-stage approach: firstly, we derive geologically relevant parameters for the size, temperature, and distribution of small intrusions in caldera settings. Specifically, we use the eroded Breiðuvík caldera in North-east Iceland and the active Krafla volcano in northern Iceland as field examples. Krafla is of particular interest due to its high potential for the exploitation of supercritical geothermal fluids [e.g. Elders et al. 2014; Scott et al. 2015]. We combine the geological observations with the results from laboratory permeability experiments on samples from Krafla by Eggertsson et al. [2020] to set up 2D finite element models that explore the influence of small, shallow intrusions on fluid flow through porous media at depths too shallow to expect super-

critical fluid. More specifically, we compare the impact of several small magmatic intrusions with that of larger intrusions. We show that networks of small intrusions emplaced at shallow depths and close in time and space may be viable targets for geothermal prospection. Our conclusions support suggestions that small intrusions may already be fuelling successful geothermal extraction sites [cf. Saubin et al. 2021].

## 2 GEOLOGICAL SETTING

We selected two field areas to extract geologically relevant observations and data on geothermal systems in caldera settings, an eroded caldera volcano (Breiðuvík) that allows insight into subsurface structures and an active caldera volcano (Krafla) with ongoing geothermal exploration and exploitation.

### 2.1 Breiðuvík volcano

The Neogene Breiðuvík volcano is part of a larger cluster of caldera volcanoes in the northern eastfjords of Iceland (Figure 1) [Burchardt 2008; Berg et al. 2014; Burchardt et al. 2022; Kennedy et al. 2022]. Large volumes of rhyolitic lavas cover basaltic lava flows that together form the volcano flanks. During an ignimbrite-forming eruption, an approximately 4 × 6 km large collapse caldera formed and was first filled with ~300–400 m of ignimbrite (Figure 1). Following the establishment of a caldera lake, the caldera depression was subsequently filled during the emplacement of subaqueous basaltic hyaloclastites and pillow breccias, before being covered by subaerial basaltic lavas. Numerous dykes and sills cut through the ignimbrite and lake sediments in the landmark Hvítserkur mountain [Kennedy et al. 2022] and the neighbouring Leirfjall. Some of the dykes connect to the overlying hyaloclastites and lavas through volcanic vents (for details see Section 4.1).

### 2.2 Krafla volcano

Krafla is an active volcano in the North Volcanic Zone of Iceland (Figure 1). It comprises an extensive fissure swarm most recently active in the historical Mývatn Fires (1724–29 CE) and Krafla Fires (1975–84 CE) eruptions, together with a 10 × 18 km collapse caldera that is mainly filled with basaltic deposits [e.g. Ármannsson et al. 1987; Mortensen et al. 2015; Kennedy et al. 2018]. While fractures associated with the fissure swarm extend 60–70 km from the caldera, eruptive fissures are confined to the nearest 7 km [Hjartardóttir et al. 2012], and the fissure swarm comprises clusters of ~3 m-wide dykes [Paquet et al. 2007]. Volcanic activity at Krafla has been predominantly bimodal, producing basaltic lava flows, hyaloclastites and tuffs, as well as rhyolitic pyroclastic deposits and domes. Geothermal exploration, syn- and post-eruptive monitoring, and scientific drilling of Krafla volcano provide an exceptionally detailed view of the shallow volcanic plumbing system that comprises some larger active magma bodies, as well as countless small intrusions [Thordarson and Larsen 2007; Kennedy et al. 2018; Árnason 2020]. The active geothermal fields at Krafla are primarily hosted within basaltic lavas and hyaloclastites, some rhyolitic ignimbrites and, at depth, gabbroic and doleritic intrusions [Ármannsson et al. 1987; Eggertsson et al. 2020]. Of critical relevance to this project is

<figure><img src="image_3.png" /></figure>
<figure><img src="image_4.png" /><figcaption>Figure 1: Geological map of the Breiðuvík caldera, East Iceland showing the location of the mountains Hvítserkur and Leirfjall. Background satellite image from Google Earth. The inset shows a map of Iceland with the location of the active volcanic zones (VZ in the N, W, and E) in dark grey shading, as well as the location of Krafla volcano and Breiðuvík caldera map. Map modified from Burchardt et al. [2022].</figcaption></figure>

that, despite the recent basaltic volcanism, there is now am-
ple evidence for shallow intrusions of rhyolite that were still
molten in 2009 [Elders et al. 2011], and derive from partial
melting of older basaltic formations [Jónasson 1994; Árnason
2020; Saubin et al. 2021].

3 METHODS

In order to produce numerical models that are informed by
observations and measurements from nature, we employ a
multi-method approach to estimate input parameters. The
geometry of the numerical models is based on field mapping,
photogrammetry, and virtual outcrop mapping. The host-

rocks for the model are chosen based on geological setup of
the Breiðuvík and Krafla calderas. Model temperatures are
derived from mineral chemistry analysis and geothermome-
try, and model rock properties are estimated from laboratory
permeability experiments performed on rocks from Krafla by
Eggertsson et al. [2020].

3.1 Geological and virtual outcrop mapping

To inform the geometric setup of our numerical models, we
mapped and sampled dykes and sills in Breiðuvík caldera in
northeast Iceland. Mapping included the on-site measurement
of intrusion orientations and thicknesses. We identified the

<figure><img src="image_5.png" /></figure>
main host-rock types, intrusion shapes, and looked for evidence of former hydrothermal circulation. We sampled representative intrusions, as well as two host-rock types (ignimbrite and lacustrine sediments). We used an Unmanned Aerial Vehicle (UAV; DJI Phantom 3 Pro) to capture aerial photographs of the steep-sided outcrops. 3D virtual outcrop models were then produced from these photographs using the software Agisoft Photoscan (now known as Metashape®) [Kennedy et al. 2022]. Ground control points were not used due to field constraints, and the resultant digital surface model is referenced using the GPS unit built into the UAV, with horizontal accuracy of ±5 m. Mapping of the outlines of the intrusion was performed using the LIME software [Buckley et al. 2019]. The mapped features of the interpreted 3D outcrops were projected onto a vertical cross-section with MOVE 2017.2™†. The relative proportion of intrusions in the section was then determined using the widely used image analysis open-source software ImageJ.

## 3.2 Mineral chemistry analysis and geothermometry

To provide an estimate of emplacement temperature of the magmatic intrusions in Breiðuvík caldera and to choose realistic modelling input temperatures, we analysed the mineral chemistry of samples from the intrusions. The analysis was performed at Uppsala University (Sweden) using the field emission source JEOL JXA-8530F Hyperprobe. The conditions for analysis were 15 kV accelerating voltage and 10 nA probe current with 10 s on peak and 5 s on lower and upper background. The detection limit for major elements is 90 ppm, and minor or trace elements, including Mn, Cr, and Ni, are not detected below 135 ppm. The beam diameter was set to 5 µm for feldspars and 1 µm for olivine and pyroxene. Uncertainties in element concentration were determined from analyses of Smithsonian Institute mineral standards. Higher element concentration gives lower uncertainties, therefore elements with concentration ≥10 wt.%, which are typically SiO₂, Al₂O₃, MgO, and CaO have uncertainties of ≤1.5 % s.d. Elements with concentrations between 5 and 10 wt.% and 2 and 5 wt.% have analytical uncertainties ≤2.2 % s.d. and ≤4.5 % s.d., respectively. Analytical uncertainties of up to 10 % s.d. apply for minor elements with concentration between 0.2 and 1.5 wt.% [Barker et al. 2015]. Pyroxene Mg# was estimated using the following equation:

$$Mg\# = \frac{Mg}{Mg + Fe_{tot}} \quad (1)$$

using the molecular proportions of Mg and Fe.

We used geothermometry on the same samples to acquire the input temperature of magmatic intrusions for the numerical modelling. Chemical analyses of elongated feldspar crystals, commonly <500 µm long and <100 µm wide, were used as input data for the plagioclase-melt equilibrium thermobarometer of Putirka [2008]. The thermobarometer requires an equilibrium test to find appropriate mineral-melt pairs. Whole-rock data of samples from the Holmar formation in eastern Iceland were used as melt composition for this study

[Óskarsson and Riishuus 2013] (Supplementary Material 1). Exchange of anorthite and albite is the deciding parameter for the equilibrium test. The mineral-melt pairs need to equal $K_D[Ab-An] = 0.28 \pm 0.11$ if $T > 1050$ °C but $K_D[Ab-An] = 0.1 \pm 0.05$ if $T < 1050$ °C to be counted in equilibrium. Mineral-melt pairs that met these requirements were then used for thermobarometric calculations. Equation 24a from Putirka [2008] was used for temperature estimates, which has standard error of estimation (SEE) of ±36 °C, and Equation 25a for pressure estimates, which has SEE of ±3.8 kbar. The H₂O content of the whole-rock basalt is needed for the thermobarometer, and 0.4 wt.% was assumed as a minimum for shallow Icelandic basalts [Nichols et al. 2002; Lehmann et al. 2007], although at depths >100 m this value can also be higher depending on the level of H₂O saturation [cf. Schopka et al. 2006].

## 3.3 Numerical Modelling

Numerical models were built using the commercially available COMSOL Multiphysics v5.6 software†. The software package uses the finite element method to solve the equations for heat transfer and fluid flow in porous media. The heat transfer module includes both conduction and convection. The fluid flow is modelled using Darcy's law with the assumption of a thermal equilibrium between the fluid and the porous medium. The host-rock around the intrusions is assumed to be homogeneous, porous, and saturated with pure water. Phase changes of the fluid (e.g. boiling) are not modelled. No recharge or outflow of water is assumed in the model domain. The bottom boundary is assigned a constant temperature, and set as impermeable with respect to fluid flow [cf. e.g. Hayba and Ingebritsen 1997]. For simplicity, intrusions are introduced as hot liquid bodies in the models with a temperature defined by our thermometry results. For details see Section 5 and Appendix D.

Every model domain was divided into up to 11,000 triangular elements. The size of elements in each model varies in three different domains depended on the geometry of the intrusions and the type of host-rock in the model. Detailed information about the meshes is provided in Appendix B. All models were run for a total of 3 × 10⁹ s (ca. 60 years), a duration that was determined empirically based on test models. The computational time steps in each model run were adjusted by the software and depend on an error tolerance, i.e. the error of calculations in each time step remains below a certain tolerance (1 in 1 million in the case of our models).

# 4 SMALL INTRUSIONS IN TWO CALDERA SETTINGS

## 4.1 Geological observations in the Breiðuvík caldera and intrusion density

Our detailed mapping and sampling of Breiðuvík caldera identified three main rock units filling the caldera that are intruded by numerous basaltic intrusions. The rock units are a >300 m-thick rhyolitic ignimbrite that is mostly unwelded in the outcrops we mapped for this study and best exposed in the mountain Hvítserkur (Figure 2A). We also identified a succession of at least 70 m of lacustrine sediments overlying the ign-

<footnote>* https://www.agisoft.com</footnote>

<footnote>† https://www.mve.com/</footnote>

<footnote>‡ https://www.comsol.com</footnote>
imbrite that had been previously mapped as hyaloclastite tuff and are best exposed in the mountain Leirfjall [Vogler 2014] (Figure 2B, D). These sediments are mostly silicic in composition and comprise sandstones, siltstones, and mudstones with diverse sedimentary features that include cross-bedding, graded bedding, ripple marks, desiccation cracks, and water escape structures, together with organic matter. On top of Hvítserkur, basaltic pillow breccias and hyaloclastites overlie a <2 m thick sedimentary succession and grade upwards into subaerial basaltic lava that forms the present-day mountain top (Figure 2A). These observations indicate that the caldera lake was mostly shallow and was eventually filled entirely by erupting basalt [Burchardt et al. 2022].

Numerous sills and dykes, as well as some agglomerate-filled vents, cut the ignimbrite and the lacustrine sediments (Figure 2). From our virtual outcrops, we calculated the proportion of intrusions to be ca. 30 % in Leirfjall (Appendix C). Sills range in thickness from <1 to 40 m, while dykes range in thickness from centimetres to ~10 m. Notably, we observed numerous examples of intrusions cross-cutting each other, without any change to the trajectory or geometry of the younger intrusion (Figure 2C). However, there are also examples of intrusions that change direction to follow alongside other intrusions (Figure 2C), forming clusters with some host-rock in between intrusions or even merging into one intrusion (e.g. Figure 2D). The thickest intrusion in the Breiðuvík caldera is a 40 m-thick sill that likely consists of four or more individual sills (Figure 2D). The different units in this sill become mappable due to differences in erosion, but in the field, no chilled margins between units can be identified.

Within both the ignimbrite and the sedimentary rocks, we found evidence of the thermal impact of the intrusions, as well as geothermal activity that was likely related to heat transfer from the intrusions (Figure 3). Close to the intrusions, the ignimbrite is often black, welded, and discoloured (Figure 3A and B) [Kennedy et al. 2022]. Notably, textural analyses of this welded ignimbrite shows a range of thermally driven textural responses with complex implications for porosity and permeability [Kennedy et al. 2022]. Nevertheless, the welded portion is more resistant to weathering than the unaffected, unwelded ignimbrite and the dykes themselves, and so the dykes are often flanked by walls of black ignimbrite (Figure 3A). Notably, these ignimbrite walls are absent for some dykes where only some centimetres of the ignimbrite are affected, but can be up to about 1 m thick for others. Discoloured ignimbrite (between centimetres and tens of metres from the intrusion) and discolouration associated with fractures extending for several metres from intrusions into the ignimbrite provide evidence for the circulation of hydrothermal fluids.

Sedimentary rocks near intrusions are metamorphosed to hornfels and appear black, densified, and more resistant to weathering than those further away from the intrusions (Figure 3C). For instance, there is a 3 m-thick zone of hornfels metamorphosed sedimentary rocks above the largest sill on Leirfjall (Figure 4A), and thermal metamorphism of the underlying 1 m of sediment, including localised melting and subsequent rheomorphic deformation. Above the sill and the hornfels-metamorphosed sediments, discoloured sedimentary

layers, mineral veins, as well as discolouration along veins, indicate alteration due to hydrothermal circulation associated with the intrusions (Figure 3D and E).

In summary, contact metamorphism and hydrothermal alteration of the host-rocks adjacent to the intrusions range in thickness from a few cm to several metres. However, establishing a correlation between intrusion thickness and the thickness of the contact metamorphic aureole will require a dedicated study.

## 4.2 Petrography, mineral chemistry, and geothermometry

### 4.2.1 Petrographic sample description

We collected samples of two representative dyke rock types from Leirfjall with different field visual characteristics, hereafter referred to as dyke A and B. Dyke A is less resistant to weathering than B and produces round weathering surfaces, while dyke B is resistant to weathering and has a blocky appearance (Figure 4A). In the outcrop, both dykes appear to crosscut the largest sill, and dyke B cuts dyke A (Figure 4A). Dyke A includes opaque minerals both as euhedral phenocrysts and in the groundmass (Figure 4B). Phenocrysts of opaque minerals are up to 1 mm in length and, together with opaque phases in the groundmass, they make up about 40–50 % of the sample. In addition to the opaque minerals, plagioclase laths up to 1 mm long constitute about 20 % of the sample. Pyroxene microlites, less than 0.1 mm in length, are evenly scattered around the sample in lower amounts than plagioclase and opaque minerals. In contrast, the groundmass of dyke B is dominated by elongated plagioclase crystals and overgrowing phenocrysts of euhedral olivine and pyroxene (Figure 4C). In addition, there are few larger (1–3 mm in diameter) phenocrysts of euhedral plagioclase in dyke B.

Samples collected from the largest sill in Leirfjall (Figure 4A) have a crystalline groundmass and a porphyritic texture (Figure 4D) similar to dyke B. Prismatic plagioclase crystals make up about 30–40 % of the groundmass. The size range of the plagioclase crystals is narrow, and the longest are about 0.5 mm. Pyroxene crystals are scattered throughout the groundmass in between the plagioclase crystals. The rest of the groundmass is composed of opaque minerals and microcrystals. Fractured olivine phenocrysts exist in most samples in low amounts and occupy up to about 5 % of some of the samples. Most of the olivine phenocrysts are of similar size and reach a diameter of approximately 1 mm.

We also collected samples from the host-rocks to the magmatic intrusions, including ignimbrite and lake sediments. The Hvítserkur ignimbrite contains pumice clasts that are several mm to <10 cm in size in a sub-millimetric matrix of glassy material (Figure 4E, F). Samples of the sedimentary host-rock collected next to magmatic intrusions are thermally altered by the intrusion and extremely fine grained with grain sizes of ≤0.1 mm (Figure 4G), implying that no or insignificant textural coarsening occurred during metamorphism. Sediments collected further away from the intrusion display larger grain sizes where the largest grains are about 0.5 mm, indicating that intrusions may have preferably invaded fine-grained host-rocks. The coarser-grained sedimentary rocks contain plagio-

<figure><img src="image_6.png" /></figure>
<figure><img src="image_7.png" /></figure>

Figure 2: Virtual outcrops produced from images taken by a UAV (drone) of the mountains [A] Hvítserkur and [B] Leirfjall with stratigraphy indicated with arrows, contacts between stratigraphic units are indicated by orange dashed lines, and yellow lines outline the intrusions. This mapping was performed using the software LIME [Buckley et al. 2019]. Dykes are mostly sub-vertical, while sills are sub-horizontal. Location of photographs in this study are indicated by white/grey boxes). [C] UAV image of dykes in ignimbrite in Hvítserkur (see location in [A]) with interpretation of the margins of individual dykes indicated by yellow lines. Note that dykes tend to either cluster or crosscut each other. [D] Drone image of intrusions in lake sediments in Leirfjall (see location in [B]) with interpretation indicated by yellow lines. Note the range of intrusion thicknesses and that the thickest sill is composed of several intrusive units.

<figure><img src="image_8.png" /></figure>
<figure><img src="image_9.png" /></figure>

Figure 3: Field photographs of observations that indicate geothermal activity associated with dykes and sills in the study area. Photograph locations are indicated in Figure 2. [A] View of dykes and a vent in ignimbrite in Hvítserkur. Note that dykes are flanked by a contact aureole of black, glassy looking ignimbrite cm to m thick, as well as discoloured ignimbrite. Red lines outline the affected ignimbrite. Dotted lines indicate the width of ignimbrite affected by the dyke in the foreground. See people for scale. The unaffected paler ignimbrite is seen in the background. [B] Example of two adjacent dykes in Hvítserkur flanked by glassy and discoloured ignimbrite. Dotted lines indicate the width of ignimbrite affected by the two marked dykes. The solid lines indicate the width of the dykes. See person for scale. The unaffected paler ignimbrite is seen in the background. [C] Hornfels of lacustrine sedimentary rock located 3 m above the thickest sill in Leirfjall. Tape measure is 0.2 m long. [D] Network of mineral veins in mudstone located within 10 m above the thickest sill in Leirfjall. Tape measure is 0.5 m long. [E] Localised veins of greenish discolouration of mudstone occur in the 10 m above the thickest sill in Leirfjall. Glove for scale.

<figure><img src="image_10.png" /></figure>
<figure><img src="image_11.png" /></figure>

Figure 4: UAV photo of intrusions in Leirfjall and photomicrographs of thin sections of intrusions and host-rocks in the study area. [A] Photograph, taken by the UAV, of the southwestern part of the summit of Leirfjall showing basaltic intrusions in layered caldera-lake sediments. White patch below dyke A is snow. Scale varies with perspective. The main features are outlined in yellow. [B] Thin section image of dyke A. Plag = plagioclase. Cross-polarised light. [C] Thin section image of dyke B. Ol = olivine, cpx = clinopyroxene. Cross-polarised light. [D] Thin section image of a sample from the thickest sill in Leirfjall. Cross-polarised light. [E] Thin section image of Hvítserkur ignimbrite displaying a pumice fragment in a matrix of glass shards. Plane-polarised light. [F] Thin section image of Hvítserkur ignimbrite displaying a mafic enclave (left middle) in a matrix of glass shards. Ves = vesicle. Plane-polarised light. [G] Thin section image of hornfels mud-siltstone sampled above the thickest sill in Leirfjall. Plane-polarised light.

<figure><img src="image_12.png" /></figure>
clase with both needle-like and tabular shapes in a matrix of indistinguishable material.

### 4.2.2 Mineral chemistry analysis

Spot analysis of samples of the magmatic intrusions was performed on four different minerals: feldspar, olivine, pyroxene, and iron-titanium oxides (see Supplementary Material 1 for complete analytical results). For an estimation of the intrusion temperature, we analysed 117 groundmass feldspar crystals with major element totals between 98 and 101 wt.%. Most of the analysed feldspars are elongate with long axes <500 µm and short axes between 10 and 100 µm. The anorthite content of the feldspars ranged from An37 to An74 with an average of An64 (Figure 5A). In a few cases, the size of feldspar crystals allowed comparison of the core and rim compositions of the crystal. In general, feldspars are normally zoned, with lower anorthite contents towards their rims.

Analysed pyroxene crystals were phenocrysts with crystal sizes >0.1 mm and in many cases overgrown by feldspar needles (see e.g. Figure 4B). The same range of acceptable analytical total was used for pyroxene crystals as for feldspar (98–101 wt.%), resulting in 39 accepted analyses. The average composition is Wo40En39Fs21, with maximum and minimum values of Wo42En47Fs29 and Wo32En29Fs15, respectively. Mg# of the pyroxenes range from 50 to 76 % with an average value of 65 %. The pyroxene crystals show narrow compositional ranges, but pyroxenes in the dykes have slightly lower Ca and higher Fe and Mg than those in the sills (Figure 5B).

The diameter of analysed olivine crystals ranges from approximately 0.1 to 0.5 mm. Most of these crystals are fractured with elongated feldspars infilling the fractures (Figure 4B, C). Olivine rims are irregular, indicating resorption.

### 4.2.3 Geothermometry

According to the $K_D[Ab-An]$ equilibrium test, 88 % (103 out of 117 data points) of the plagioclase crystals were in equilibrium with the chosen whole-rock data (see Appendix A) [Öskarsson and Riishuus 2013] and therefore used for the geothermometry (Figure 6). Temperature values from the thermometry calculations ranged from 1120 to 1139 °C with an average value of 1130 °C (Figure 6). We therefore used 1130 °C as the initial temperature of the magmatic intrusion in the numerical models. Pressures were not used as model input parameters because the error associated with the geobarometers is high (SEE = ±3.8 kbar). Instead, the intrusion depths in the models (ca. 1–1.5 km) are based on geological constraints from the Krafla caldera, i.e. small sheet intrusions are abundant at shallow depths [Scott et al. 2015; Kennedy et al. 2018].

## 5 MODEL SETUP BASED ON CONSTRAINTS FROM THE GEOLOGY OF BREIÐUVÍK AND KRAFLA CALDERAS

As this study aims to use geologically informed numerical models, this section will explain how we set up the models based on constraints we derive from our field examples.

### 5.1 Model geometry

The geometry of our numerical models is informed by our field observations from Breiðuvík and previous scientific and geothermal exploration at Krafla caldera [cf. Weisenberger et al. 2015; Árnason 2020]. We chose the model dimensions and intrusion depths according to borehole depths into geothermal systems in Krafla [e.g. Árnason 2020], while intrusion geometries and dimensions are based on observations in Breiðuvík, albeit strongly simplified to account for a feasible model resolution and to test our research question regarding the geothermal potential of small intrusions. The models are 2D, 2500 m-wide and 1500 m-high/deep to avoid boundary effects (Figure 7). Intrusions were placed on the left side of the model. This left side is defined as a vertical symmetry plane, which allows us to save computational resources by modelling only half of the domain.

The shape of intrusions is highly simplified compared to nature, as host-rock contacts are assumed to be planar and intrusion ends rounded. The thickness of the intrusions is either 40 m or 5 m, based on the thickness of major intrusions observed in the Breiðuvík caldera (see Section 4.1). To test whether one large intrusion has the same effect on the geothermal circulation as several smaller ones, we designed models with either a single 40 m-thick sill/dyke or eight 5 m-thick sills/dykes with 10 m spacing (Figure 7), i.e. the cumulative thickness of the thin intrusions is equal to the thickness of the thick intrusion. These intrusion arrays represent the clusters of parallel intrusions observed in the Breiðuvík caldera. The length of the intrusions in the models is kept constant at 500 m to ensure comparability. Hence, the four different intrusion arrays modelled are: a) a single, thick sill, b) eight thin sills, c) a single, thick dyke, and d) eight thin dykes (Figure 7).

### 5.2 Model temperatures

The initial temperature of the host-rock was set to follow a vertical temperature gradient of 200 °C km⁻¹ in the top 1 km to simulate a high-temperature (i.e. high-enthalpy) geothermal system. Below a depth of 1 km, the initial temperature of the host-rock was assumed to be 200 °C, because temperature profiles from high temperature areas show a tendency towards more stable temperatures around 1 km depth [Arnórsson et al. 2008]. The initial pressure in the system was defined as hydrostatic. The surface temperature at the top boundary was set to 0 °C and pressure to 1 atm (101,325 Pa). In order to achieve a stable background temperature distribution, including the development of hydrothermal convection in the host-rock, we ran the models for 10⁹ s (31.7 years) before adding the intrusion(s).

The intrusions are introduced instantaneously with an initially uniform temperature of $T = 1130$ °C. We neglect thermal convection in the intrusions as its extent and significance is still debated [Marsh 1989; Eichelberger 2020]. Given the nearly isothermal conditions and the small spatial extent of the convective layer, convective heat transport is generally low, and purely conductive heat transport provides a good estimate of the total heat transport (as evidenced by a Nusselt number around unity; see Figure 9 in Marsh [1989]). For details see Appendix D.

<figure><img src="image_13.png" /></figure>
<figure><img src="image_14.png" /><figcaption>Figure 5: Results from mineral chemistry analysis of dykes (red circles and orange squares) and sills (blue triangles) collected from Leirfjall. [A] Compositional classification of feldspars. [B] Compositional classification of pyroxenes.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 6: Results from equilibrium test and thermobarometry for plagioclase mineral in three representative intrusions (dykes and sills) in Leirfjall.</figcaption></figure>

<table><caption>Table 1: Permeability of the magmatic intrusions that changes linearly as a function of temperature [cf. Hanson and Barton 1989; Foumier 1991; Hayba and Ingebritsen 1997; Gerdes et al. 1998].</caption><thead><tr><td>Temperature T [°C]</td><td>Permeability k [m²]</td></tr></thead><tbody><tr><td>360</td><td>1.55 × 10<sup>-16</sup></td></tr><tr><td>500</td><td>1 × 10<sup>-17</sup></td></tr><tr><td>750</td><td>5 × 10<sup>-18</sup></td></tr><tr><td>900</td><td>1 × 10<sup>-22</sup></td></tr></tbody></table>

Since we consider the intrusions to be impermeable at temperatures above the solidus, we set permeability to $1 \times 10^{-22} m^2$ at $T > 900 °C$ [e.g. Gerdes et al. 1998]. The permeability then increases linearly with decreasing temperature to account for e.g. fracture formation in the intrusion (see Table 1) [cf. Foumier 1991; Hayba and Ingebritsen 1997].

## 5.3 Host-rock properties

The rocks hosting the intrusions in our study areas include basaltic lavas, hyaloclastites, ignimbrites, and lacustrine sediments (see Section 4.1) [Mortensen et al. 2014; Eggertsson et al. 2020]. Each of these rock types is known to have a wide range of physical and mechanical properties [e.g. Schaefer et al. 2015; Eggertsson et al. 2020; Heap et al. 2020a; Weaver et al. 2020; Di Muro et al. 2021; Heap and Violay 2021], due to composition, eruption type and style, and alteration, amongst other factors [as discussed by Heap and Violay 2021]. It is also of note that the response of volcanic rocks to thermal stimulation is complex and temperature-dependent [Mordensky et al. 2019; Kennedy et al. 2022]. Similarly, lithologically distinct host-rocks are known to respond differently to thermal stressing from intrusions [Saubin et al. 2019]. The thermal properties (the thermal conductivity and heat capacity) of our model host-rocks are listed in Table 2 and considered to be constant for the following reasons: preliminary modelling in preparation of this study revealed that permeability is the principal parameter controlling fluid flow in the host-rock surrounding small intrusions, and the fact that the effect of the thermal properties is minor to negligible in our models in contrast to the complex thermo-chemical changes seen in the field and in experiments [Saratovich et al. 2015; Mordensky et al. 2019; Eggertsson et al. 2020].

Based on the results of these preliminary, exploratory models, we therefore chose to use two simplified host-rock types only (Table 2): one with low permeability and one with high permeability [cf. Árnason 2020]. The low-permeability host-

<figure><img src="image_16.png" /></figure>
<figure><img src="image_17.png" /><figcaption>Figure 7: Sketches of the model geometries with intrusions in red and host-rock in grey. Yellow circles mark the locations of temperature measurements in Figure 10. Note that the drawings are not exactly to scale to allow visibility of small features. 1 atm = 101,325 Pa.</figcaption></figure>

rock represents rocks such as dense basaltic lava, welded ignimbrite, and mudstone, while the high-permeability host-rock represents rocks such as vesicular lavas, hyaloclastite, and sandstone. Permeability values for the host-rocks are based on laboratory experiments performed on rocks from Krafla volcano by Eggertsson et al. [2020]. To account for the effect of overburden pressure on host-rock permeability, we use the empirical results of permeability as a function of effective pressure by Eggertsson et al. [2020] as a starting point. The permeability variations with effective pressure observed by Eggertsson et al. [2020] can be fitted to an exponential function, which we converted to variations with depth assuming lithostatic conditions for the host-rock and hydrostatic conditions

for the pore fluids (Figure 8 and Appendix D). The resulting host-rock permeabilities are similar to the ones assumed by Árnason [2020] (high-permeability rock with a permeability of $4.9 \times 10^{-14} \text{ m}^2$, and low-permeability rock with a permeability of $9.8 \times 10^{-16} \text{ m}^2$).

As our model host-rocks are based on laboratory experiments, the two basic rock models for permeability cannot account for fractures larger than the scale of the laboratory samples (>mm to 2–3 cm), which can considerably increase the permeability of a rock mass [e.g. Ebigbo et al. 2016; Hobé et al. 2018; Saubin et al. 2019] (see Section 6.2). In order to investigate the possible effect of fractures, we introduce a third host-rock type with elevated permeability (“super-high permeability”).

<figure><img src="image_18.png" /></figure>
ability"). As long as there are several interconnected fractures
and the flow through the fractures is laminar, the Darcy flow
model of porous flow is applicable. Of course, it is possi-
ble that the geothermal flow is strongly controlled by a sin-
gle, large fracture, which would lead to strongly localised and
anisotropic flow. However, as we do not have any such data
about the host-rock in the study area and as we did not ob-
serve any large veins, we deem it most appropriate to model
fractured rock units in terms of isotropic Darcy flow with in-
creased permeability.

The combination of the four intrusion geometries and three
host-rock permeabilities results in a total of 12 numerical mod-
els, as listed in Table 3.

<figure><img src="image_19.png" /><figcaption>Figure 8: Permeability as a function of depth for the three model host-rocks, one with high, one with low, and one with super-high permeability. The functions are derived from fitting exponential functions to the empirical values derived by Egertsson et al. [2020]. See also Appendix D, Figure D3, D4, D5, D6, and D7</figcaption></figure>

Our modelling results show the cooling of all intrusions,
as well as the effect of the intrusions on the temperature field
in the surrounding rock (Figure 9). Regarding the temperature
effect on the host-rock, we observe that conductive heat trans-
port is dominant in all models with low-permeability host-
rocks, while heat from the intrusions in the models with high
and super-high permeability host-rocks triggers vigorous con-
vection above the intrusions.

Cooling of the intrusions at their tips is on the order of
decades (Figure 10). While all intrusions experience a rapid
loss of heat by several hundred degrees immediately after
emplacement, cooling considerably slows down within one
to five years. The single large intrusions cool on longer
timescales than the arrays of small intrusions, with the single
dyke cooling slowest in both the low- and high-permeability
models. The single dyke low permeability (1D-LOW) model
displays a temperature stabilisation between 700 and 600 °C

for the three decades after emplacement, a similar rapid drop
and stabilisation is seen with the array of eight dykes (8D) and
eight sills (8S) in low permeability host-rocks stabilising at a
temperature between 450 and 500 °C and between 350 and
400 °C, respectively (Figure 10). The single large sill displays
a more continuous heat loss and stabilises only after more
than 20 years at temperatures above 400 °C. The arrays of
small intrusions experience a considerably higher initial heat
loss than the large intrusions and their temperatures stabilise
at about 50 to 150 °C below the large intrusions for all model
host-rocks. Generally, the lower the host-rock permeability,
the slower the intrusions cool. However, three decades after
emplacement, all intrusions are still considerably above the
background temperature (Figure 10).

The temperature field of the models highlights that all in-
trusions create a narrow thermal aureole by heating their im-
mediate surroundings within the modelled time (Figures 9
and 11). Here, we define the thermal aureole as the area in the
host-rock in which temperatures are increased with respect
to the background geotherm due to the presence of the intru-
sions, but excluding the plumes of hot fluid rising above the
intrusions in some models. The width of the thermal aureole
is primarily controlled by the host-rock permeability, with the
widest aureoles for the low permeability host-rock (Figure 11).
In those models, the thermal aureole is approximately 12 % of
the sill width (60 m) for the single sill, and increasingly wider
for the sill array (ca. 80 m), single dyke (ca. 100 m), and dyke
array (ca. 110 m). The relationship of the aureole widths is the
same for the high and super-high permeability host-rocks (1S
< 8S < 1D < 8D), but absolute widths are much smaller, i.e.
approximately 5 to 50 m for the high permeability models and
1 to 25 m for the super-high permeability models. Notably,
in the models with super-high permeability, temperatures in
the thermal aureole drop beneath the background temperature
because of the downwelling of cold water (Figure 11).

Our models also highlight that the intrusion geometry has
some influence on how intrusions cool. While dykes cool from
their upper tip downwards, sills cool from their lateral ends
inwards. Heating of the host-rock in between intrusions in the
arrays creates a common thermal areole for the entire array
within the first years after intrusion emplacement. Once the
rock between the thin intrusions is heated, the array acts as
one heat source and the rapid cooling of the array slows down
considerably (cf. Figure 10). Interestingly, the intrusion arrays
create wider thermal aureoles than the single, thick intrusions
in all models (Figure 11). The observation that dykes create
wider thermal aureoles than sills may be explained by the
position of the thermal plumes above the intrusions and the
prolonged upward transport of hot water along the dykes from
greater depths (cf. Figure 9).

6 DISCUSSION

6.1 Discussion of geological observations

Our observations from Breiðuvík caldera and previously doc-
umented geological information from Krafla stress the com-
plexity of caldera volcanoes as hosts of geothermal reservoirs.
Not only do calderas contain a multitude of different host-

<figure><img src="image_20.png" /></figure>
<table><caption>Table 2: Overview of rock physical properties used in the modelling performed for this study. In order to simplify the model input parameters, we use three rock types only, i.e. a low-permeability, a high-permeability, and a super-high permeability host-rock. Densities are based on a bulk dry density that assumes a solid matrix density of 2900 kg m⁻³. Porosities are derived from the range of rocks measured by Eggertsson et al. [2020], thermal conductivity is from Robertson and Peck [1974] and based on dry basalt; heat capacity from Heap et al. [2020b].</caption><thead><tr><th>Host-rock (HR) type</th><th>Solid-rock density [kg m<sup>-3</sup>]</th><th>Porosity</th><th>Permeability at the Earth's surface [m<sup>2</sup>]</th><th>Thermal Conductivity [W m<sup>-1</sup> °C]</th><th>Heat Capacity [J kg<sup>-1</sup> °C]</th></tr></thead><tbody><tr><td>Low permeability HR</td><td>2850</td><td>0.03</td><td>1.32 × 10<sup>-15</sup></td><td>1.75</td><td>1050</td></tr><tr><td>High permeability HR</td><td>2850</td><td>0.28</td><td>7.31 × 10<sup>-14</sup></td><td>1.75</td><td>1050</td></tr><tr><td>Super-high permeability HR</td><td>2850</td><td>0.35</td><td>10<sup>-12</sup></td><td>1.75</td><td>1050</td></tr></tbody></table>

<figure><img src="image_21.png" /><figcaption>Figure 9: The modelled temperature field ca. 30 years after intrusion. Note that a 600 m-wide section of the models is shown only. The colour scale focusses on the temperature range interesting for geothermal exploration. Temperatures above 300 °C do occur. See Table 3 for a list of all models. See text for description.</figcaption></figure>

rocks at various stages of alteration and fracturing, they are also intruded by a complex network of variably sized intrusions. Our geological investigations of Breiðuvík caldera provide qualitative and quantitative constraints on dyke and sill networks cutting through diverse host-rocks. Similar relation-

ships for caldera intrusions are well documented by geothermal and scientific drilling at Krafla caldera [e.g. Mortensen et al. 2015; Weisenberger et al. 2015] and elsewhere in Iceland [e.g. Saubin et al. 2019; Liotta et al. 2020]. Equivalent constraints on intrusion geometry and spacing derive from stud-

<figure><img src="image_22.png" /></figure>
<figure><img src="image_23.png" /><figcaption>Figure 10: Cooling of the intrusions in all models, measured 20 m from the tip of the (top) sill and (middle) dyke (see locations in Figure 7). Measurement point was chosen for comparability. Colour refers to intrusion geometry (where blue and red curves refer to dykes and sills, respectively), line style refers to host-rock permeability. See Table 3 for explanation of the model names. Time is given relative to the onset of intrusion, in years. The fluctuations in temperature for models 8D-SHI and 8S-SHI between 25 and 30 years is due to the downwelling of cold plumes in the geothermal system.</figcaption></figure>

<table><caption>Table 3: List of the numerical models performed for this study, indicating their intrusion geometry (dyke or sill), number of intrusions (one or eight), and the host-rock (high, low, or super-high permeability).</caption><thead><tr><th>Model ID</th><th>Geometry</th><th>Host-rock</th></tr></thead><tbody><tr><td>1S-Low</td><td rowspan="3">Single thick sill</td><td>Low-permeability</td></tr><tr><td>1S-Hi</td><td>High permeability</td></tr><tr><td>1S-SHi</td><td>Super-high permeability</td></tr><tr><td>8S-Low</td><td rowspan="3">Eight thin sills</td><td>Low-permeability</td></tr><tr><td>8S-Hi</td><td>High permeability rock</td></tr><tr><td>8S-SHi</td><td>Super-high permeability</td></tr><tr><td>1D-Low</td><td rowspan="3">Single thick dyke</td><td>Low-permeability</td></tr><tr><td>1D-Hi</td><td>High permeability</td></tr><tr><td>1D-SHi</td><td>Super-high permeability</td></tr><tr><td>8D-Low</td><td rowspan="3">Eight thin dykes</td><td>Low-permeability</td></tr><tr><td>8D-Hi</td><td>High permeability</td></tr><tr><td>8d-SHi</td><td>Super-high permeability</td></tr></tbody></table>

ies of dykes and fractures associated with fissure swarms in the active and eroded rift zones of Iceland [e.g. Paquet et al. 2007; Hjartardóttir et al. 2012] and other eroded volcanoes in Iceland and worldwide [e.g. Walker 1999; Khoei et al. 2015; Karaoğlu et al. 2016; Burchardt et al. 2018; Mordensky et al. 2018; Bazargan and Gudmundsson 2019; 2020; Ruz et al. 2020; Satow et al. 2021].

Notably, our field investigations highlight examples of cross-cutting intrusions that do not change trajectory and the clustering of individual intrusions (see Section 4.1). The former observation is generally accepted as evidence of complete solidification of the earlier emplaced intrusion, while the latter is interpreted as evidence for the interaction of intrusions emplaced in a rapid manner [cf. Burchardt 2008; Chanceaux and Menand 2014]. Clustering of intrusions in the same dyke swarm has also been documented by Paquet et al. [2007] and attributed to weakening of the host-rock by hydrothermal alteration or dyke emplacement. In Iceland, rifting episodes, such as the Mývatn and Krafla fires in Krafla [Einarsson 1991; Buck et al. 2006; Sturkell et al. 2008] and the currently ongoing activity on the Reykjanes peninsula highlight that several sheet intrusions can be emplaced in the same volcanic system within weeks to years [Cubuk-Sabuncu et al. 2021; Hobé et al. 2021].

<figure><img src="image_24.png" /></figure>
<figure><img src="image_25.png" /><figcaption>Figure 11: Cooling of the intrusions and heating of the host-rocks along temperature profiles through the models at 1020 m depth during different time steps after intrusion, in years. [A]–[D] Temperature along the profile at different time steps. [E]–[F] The lateral extend of the intrusion(s) and the maximum extend of the area in which temperatures are increased because of the intrusion(s) (i.e. the thermal aureole).</figcaption></figure>

Our geological observations also highlight that small magmatic intrusions can cause contact metamorphism and circulation of geothermal fluids in fracture networks in their immediate surroundings (e.g. Figure 3). The thickness of thermal aureoles we observed in the field were similar to, or thinner than, the intrusion thickness. The lack of extensive melting

of the host-rocks in the aureole indicates that temperatures were lower than those achieved in the numerical models. This seeming mismatch can be explained by the shallow depth of our field area (<400 m). In the models, the temperature field is a sum of the intrusion heat and the background geothermal gradient at depths of up to 1500 m. Therefore, the tempera-

<figure><img src="image_26.png" /></figure>
tures in the models are not temperatures we expect to have occurred in the field examples.

Quantitative mapping of intrusion networks like the one in Breiðuvík caldera highlight the complexity of such geological systems regarding intrusion geometry, dimensions, interaction with each other, pre-existing structures, and the host-rock(s). While this complexity will necessitate further investigation to explore the geothermal viability range, the modelling results indicate that geothermal production would already be viable for these baseline cases. By informing and validating our numerical models with geological information, these baseline cases will thus reduce exploration risks and open up possibilities for exploiting shallow intrusions

## 6.2 Model limitations and potential for future studies

Although the numerical models presented in this study are partially constrained by field observations of intrusion geometries and dimensions, laboratory permeability measurements [Eggertsson et al. 2020], and geothermometry estimates from intrusions as presented in this study, the models remain highly simplified approximations of natural systems. Permeability can range over several orders of magnitude in a single rock type such as ignimbrite, especially in hydrothermal environment [Heap et al. 2019; Kennedy et al. 2022]. The model geometry of intrusions, as well as the homogeneity of the host-rocks, are oversimplified. As evident from our own observations in the field, intrusion geometries and networks are much more complex, e.g. including dykes, sills, and other intrusion geometries with different thickness and different spacing in the same outcrop. Moreover, intrusions have complex three-dimensional shapes that would yield more complex cooling patterns when modelled in 3D. Our models provide a first approximation of the thermal impact of small intrusions, and moving towards more complex, 3D scenarios requires careful sensitivity studies of all involved parameters. Future 3D modelling can then also take into consideration different orientations of intrusions and other features, as indeed, fractures within Krafla caldera are aligned along a number of different orientations that are not parallel with the main fissure swarm [Hjartardóttir et al. 2012]. This contrasts with the more uniform orientation and dimensions of fractures and dykes associated with Krafla and other fissure swarms [e.g. Paquet et al. 2007; Hjartardóttir et al. 2012; Ruz et al. 2020]. Nonetheless, the predominant fissure-parallel N to NNE orientation of fractures opening within the Krafla caldera in the 1975–1984 rifting episode indicates that the local caldera-related stress field does not always strongly influence dyke orientations [Hjartardóttir et al. 2012], and so temporal changes in caldera-related stress fields must also be an important factor.

We modelled the instantaneous emplacement of intrusions even in the arrays, corresponding to cases where successive intrusions have insufficient time to completely cool before the next intrusion is formed. Our field observations (cross-cutting relationships) and observations of rifting episodes, such as the nine-year-long Krafla Fires [e.g. Einarsson 1991; Buck et al. 2006; Sturkell et al. 2008], provide indirect and direct evidence for emplacement of sheet intrusions within the same array over periods of months to years. A successive emplacement of

several intrusions should have a strong effect on the nature and longevity of temperature anomalies in their surroundings [cf. Annen 2011; Annen 2017], as well as influencing the local stress field [Paquet et al. 2007]. Further systematic analysis of the spacing of fossil dissected intrusions [e.g. Paquet et al. 2007; Ruz et al. 2020] and the timing of intrusion events during instrumentally observed rifting episodes [e.g. Björnsson et al. 1979; Ágústsdóttir et al. 2016] will provide useful new constraints.

The use of an exponential function fitted to laboratory permeability data for Krafla rocks with low and high permeability by Eggertsson et al. [2020] allowed us to estimate the permeability of rocks at depth (Figure 8). The values from the Eggertsson et al. [2020] study cover intact rocks, as well as rocks with fractures on the scale of laboratory investigations (mm to a few cm). The range of permeabilities we chose to derive the permeability as a function of effective pressure therefore covers sample-scale fractures. However, large-scale fractures are known to significantly increase the equivalent permeability of a volcanic rock mass [e.g. Heap and Kennedy 2016; Hobé et al. 2018]. The inclusion of fractures in the models is problematic as no representative elementary volume exists for fractured rock. We therefore explored the influence of an increased permeability due to fractures by introducing a model host-rock with super-high permeability (one order of magnitude higher than the high permeability host-rock). Nevertheless, fracture networks need further investigation as they provide preferential flow-paths, and can produce large ranges of permeability values due to multi-pathing, especially near the percolation threshold [e.g. Ebigbo et al. 2016; Hobé et al. 2018, and references therein]. The potential healing or sealing of fractures represents an additional complication. We therefore neglect the influence of fractures beyond the simple permeability increase observed in Eggertsson et al. [2020]. Because of this, our results can both be based on underestimated or overestimated permeabilities of our host-rocks, which can result in an under- or overestimation of the modelled cooling timescales.

As a result of our approach, the effect of host-rock permeability decreases with depth in our models. However, we assume constant values for the thermal properties of the model host-rocks. Since porosity likely decreases as a function of depth [e.g. Wyering et al. 2014], it is likely that, for example, the thermal conductivity of volcanic rocks will also increase as a function of increasing depth [see, e.g. Heap et al. 2020a]. With more data, future modelling could also account for the depth dependence of thermal properties.

Furthermore, our geological observations and other studies highlight the influence of host-rock layering in the emplacement of intrusion [e.g. Mathieu et al. 2015]. As discussed above, it is well established that lithological changes can introduce variations in host-rock properties that can influence fluid flow and associated heat transport, which is why future modelling also needs to address these parameters.

Considering that our models are strongly simplified and that there are near limitless possibilities to expand the investigated parameter space, we would like to stress that our models simply aim to investigate the effect of the intrusions as heat sources and to compare larger, single intrusions with

<figure><img src="image_27.png" /></figure>
arrays of smaller intrusions. In nature, geothermal reservoirs would require further components, such as reservoir and cap rocks. For instance, separation of fluid circulation into shallow and deep thermal regimes at Vítismór-Leirbotnar within the Krafla geothermal system has been attributed to either a clay-rich cap rock or depth-dependent permeability [Arnason 2020], and more refined modelling of formation permeability with depth and position is required.

## 6.3 Implications for the geothermal potential of small intrusions in active caldera volcanoes: small intrusions, big effects

Our geologically informed models highlight that even intrusions much smaller than the normally considered plutons cause a locally significant increase of the sub-surface temperature, and this temperature increase occurs within decades following emplacement. Still, this timescale remains highly relevant to Krafla (where the last magmatic intrusions spanned 1975–1984 with still active magma bodies in 2009) and other frequently-active volcanic-geothermal systems. In fact, the high temperatures encountered by boreholes in Krafla indicate a potential proximity to recent sheet intrusions [Saubin et al. 2021].

Not surprisingly, our model results indicate that larger intrusions provide more heat to their surroundings for a longer time. However, the limited extent of the thermal aureole makes them a difficult target for drilling, in particular in the case of the dykes that have low thicknesses. The models also show that despite the faster cooling of arrays of small intrusions, such networks, and their thermal aureole, can act as significant heat sources for several decades following their emplacement. In fact, the magmatic-hydrothermal interactions triggered by the emplacement of several sheet intrusions may explain thermal anomalies measured several years prior to eruptions (thermal unrest) [Girona et al. 2021].

As highlighted by our field examples, the Breiðuvík and Krafla calderas, such networks of small intrusions are abundant at shallow depths. Our models also highlight that several small intrusions forming close to each other in time and space may afterwards act as a single heat source once the rock between intrusions is heated, as the intrusion array distributes the heat over a larger area, and the cooling of the outer intrusions buffers the cooling of the inner intrusions. The thermal anomaly created by such an array of near-contemporaneous small intrusions could be larger in area than that of a single thicker intrusion, and the larger thermal aureole around the arrays implies that higher volumes of water are heated by these arrays compared to the single intrusions. Moreover, in nature, the presence of both sills and dykes in the same array likely increases the thermal anomaly compared to the scenarios modelled in this study. Hence, a network of dykes and/or sills emplaced during a single rifting event may provide sufficient heat to be economically exploitable for electricity production from water heated >250 °C. While all modelled scenarios would provide sufficient heat for energy production, the intrusion arrays would be more efficient at producing high volumes of water at optimal temperatures for electricity production.

We suggest costs and risks of geothermal projects could be reduced in many aspects by exploiting networks of sheet intrusions at shallow depths. Higher resolution exploration methods should be used to improve resource assessment and drill targeting and thus reducing risk. Boreholes to shallow intrusions would be much cheaper than deep boreholes as drilling costs grow exponentially with depth. Increased flow-rates can be achieved due to higher permeabilities at shallow depths, thus allowing faster heat extraction. A power plant could thus target multiple intrusions individually. Wells that stop producing high enough temperatures for direct electricity production (ca. <150 °C) could still be used in multiple ways:

1. Super-heating the fluids using another source [Borsukiewicz-Gozdur 2010] would still allow electricity production (e.g. mixing of >300 °C fluids with 150 °C fluids or through burning biomass to get temperatures relevant for electricity production);

2. Cascade heating [Rubio-Maya et al. 2015] would allow direct heat usage until temperatures reduce to ~18 °C;

3. Reinjection of geothermal fluids produced elsewhere could allow extraction without the need of an additional injection well;

4. Energy storage (heat and other sources); and

5. CO₂ sequestration.

A system of multiple wells connected to shallow heat sources at different temperatures thus reduces risk as it allows for flexibility. This is especially true if added to an existing system based on a larger heat source, as most of the required systems would already be in place.

Our results provide vital information for future geothermal exploration or research-driven magma drilling endeavours such as the Krafla Magma Testbed (KMT) [Eichelberger 2019]. Prior to significant investment in infrastructure, such projects will require reliable geophysical identification of magma body locations and dimensions, in order to constrain the intruded volume and nature of the host-rock, and thus estimate the lifetime of potentially useful heat output. Small intrusions are admittedly hard to identify with the current geophysical methods [cf. Rooyakkers et al. 2021]. However, in the case of dyke intrusions identified by real time monitoring, additional constraints on dyke widths, which are a key influence on thermal longevity, may derive from geodetic data [e.g. Sigmundsson et al. 2014] and passive seismic data [Cubuk-Sabuncu et al. 2021; Hobé et al. 2021]. In case of sills or sill arrays, horizontal drilling may be preferred, as the geometry of sills implies that they may act as permeability barrier. Once the challenges of intrusion identification have been overcome, abundant small intrusions at shallow depths in active caldera volcanoes could become a worthwhile target for future geothermal exploration. In particular, small intrusions emplaced into low-permeability host-rocks, such as dense basalt lava, cool slowly because heat transfer to the host-rock is dominated by conduction. Natural permeability enhancement from the intrusion [Mordensky et al. 2018] or induced via mechanical hydrofracturing [e.g. Siratovich et al. 2015] or chemical stimulation [e.g. Farquharson

<figure><img src="image_28.png" /></figure>
et al. 2019] could herald efficient heat extraction from small intrusions in such host-rocks.

# 7 CONCLUSIONS

We analysed the geothermal potential of small magmatic sheet intrusions in caldera settings by studying the eroded Breiðuvík caldera in Eastern Iceland and its exposed network of solidified sheet intrusions as an analogue for the active Krafla volcano, Northeastern Iceland. Sheet intrusions of various sizes (up to 40 m thick) and geometries (dykes and sills) occur in networks that make up about 30 % of the outcrop. Field evidence suggests some of the intrusion arrays in the outcrops were emplaced over short time spans (months to years), and the host-rocks show signs of former geothermal activity.

Together with geothermometric estimates of the intrusion temperature in Breiðuvík sills and dykes, we used these measurements and our observations to model cooling and heat transport in 2D numerical models that simulated fluid flow in porous media. With respect to the geothermal potential of such intrusions, the results of these models indicate that although single, large intrusions have a stronger thermal impact on the surrounding rocks, networks of contemporaneous small intrusions may heat up a larger area in the subsurface, which implies larger volumes of high-temperature fluids. All models show that the intrusion-related thermal anomalies remain for decades, for all model host-rocks. Large volumes of >250 °C fluids may thus be available for extraction for several decades after the emplacements of intrusions. Given that borehole costs increase and geophysical imaging resolution decrease with depth, and that networks of small intrusions occur abundantly at shallow depth in active caldera volcanoes, shallow intrusions should be considered as a target for geothermal exploration.

# AUTHOR CONTRIBUTIONS

S.B., C.H., M.B. and E.B.G. formulated the research question. Numerical modelling was performed by M.B., with inputs and help from E.R., A.H., C.H., E.BG. and S.B. E.B.G. was awarded funding for fieldwork for this study under supervision of S.B. by Landsvirkjun. Field work, sampling, and photogrammetry were carried out by E.B.G., S.B., H.T., B.K., J.D., and E.S. Sample preparation and analysis, as well as petrography, mineral chemistry and thermometry analyses were performed by E.B.G. Mapping of intrusions on virtual outcrops was performed by S.B. and E.B.G. M.H. provided input on the physical properties of rocks and helped validate the model setup. A.H. contributed with advice on the geothermal potential and economic aspects, as well as the model limitations. The manuscript was assembled by S.B. and E.B.G. with input from all authors.

# ACKNOWLEDGEMENTS

The authors acknowledge a research stipend by Landsvirkjun to E.B.G. that partly funded fieldwork in Breiðuvík and advice on mineral chemistry and thermometry by Abigail Barker. S.B. was supported by a Knut and Alice Wallenberg fellowship. H.T. was supported by Royal Society University Research Fel-

lowship UF140716; H.T. and M.H. were supported by Royal Society International Exchanges grant IE150771. E.R. thanks the Swedish National Infrastructure for Computing (SNIC) at UPPMAX, partially funded by the Swedish Research Council through grant agreement no. 2018-05973, for the resources in project SNIC 018/8-219 that supported the preliminary phase of the models. M.H. was partly supported by the Interdisciplinary Thematic Institute (ITI) "Geosciences for the Energy System Transition" of the University of Strasbourg, funded under the Excellence Initiative Program. M.H. also acknowledges a Partenariats Hubert Curien (PHC) ULYSSES grant. B.K. acknowledges MBIE Catalyst grant Energy from Magma. The authors thank the editor John Browning and two anonymous reviewers for their constructive reviews of the manuscript.

# DATA AVAILABILITY

All relevant data is provided in the appendices (Appendix A, B, C, and D) and Supplementary Material alongside this article. This includes analytical results (Supplementary Material 1) and animations of temperature field simulations (Supplementary Material 2). The module recipe used with this COMSOL model, as well as simulations of the model temperature fields are freely available on the YouTube channel (https://www.youtube.com/channel/UCkcP0IZ09vTlp1FfIPg_8rg/about).

# COPYRIGHT NOTICE

© The Author(s) 2022. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

# REFERENCES

* Ágústsdóttir, T., J. Woods, T. Greenfield, R. G. Green, R. S. White, T. Winder, B. Brandsdóttir, S. Steinthórsson, and H. Soosalu (2016). "Strike-slip faulting during the 2014 Bárðarbunga-Holuhraun dike intrusion, central Iceland". *Geophysical Research Letters* 43(4), pages 1495–1503. ISSN: 1944-8007. DOI: 10.1002/2015g1067423.

* Annen, C. (2011). "Implications of incremental emplacement of magma bodies for magma differentiation, thermal aureole dimensions and plutonism-volcanism relationships". *Tectonophysics* 500(1-4), pages 3–10. DOI: 10.1016/j.tecto.2009.04.010.

* Annen, C. (2017). "Factors Affecting the Thickness of Thermal Aureoles". *Frontiers in Earth Science* 5. ISSN: 2296-6463. DOI: 10.3389/feart.2017.00082.

* Ármannsson, H., A. Gudmundsson, and B. S. Steingrímsson (1987). "Exploration and development of the Krafla geothermal area". *Jökull* 37, pages 13–30.

* Árnason, K. (2020). "New Conceptual Model for the Magma-Hydrothermal-Tectonic System of Krafla, NE Iceland".

<figure><img src="image_29.png" /></figure>
Geosciences 10(1), page 34. ISSN: 2076-3263. DOI: 10.3390/geosciences10010034.

Arnórsson, S., G. Axelsson, and K. Sæmundsson (2008). "Geothermal systems in Iceland". *Jökull* 58, pages 269–302.

Barker, A. K., V. R. Troll, J. C. Carracedo, and P. A. Nicholls (2015). "The magma plumbing system for the 1971 Teneguía eruption on La Palma, Canary Islands". *Contributions to Mineralogy and Petrology* 170(5-6). ISSN: 1432-0967. DOI: 10.1007/s00410-015-1207-7.

Bazargan, M. and A. Gudmundsson (2019). "Dike-induced stresses and displacements in layered volcanic zones". *Journal of Volcanology and Geothermal Research* 384, pages 189–205. ISSN: 0377-0273. DOI: 10.1016/j.jvolgeores.2019.07.010.

- (2020). "Stresses and displacements in layered rocks induced by inclined (cone) sheets". *Journal of Volcanology and Geothermal Research* 401, page 106965. ISSN: 0377-0273. DOI: 10.1016/j.jvolgeores.2020.106965.

Berg, S. E., V. R. Troll, S. Burchardt, M. S. Riishuus, M. Krumbholz, and L. E. Gústafsson (2014). "Iceland's best kept secret". *Geology Today* 30(2), pages 54–60. ISSN: 0266-6979. DOI: 10.1111/gto.12042.

Björnsson, A., G. Johnsen, S. Sigurdsson, G. Thorbergsson, and E. Tryggvason (1979). "Rifting of the plate boundary in north Iceland 1975–1978". *Journal of Geophysical Research* 84(B6), page 3029. ISSN: 0148-0227. DOI: 10.1029/jb084ib06p03029.

Björnsson, G. (1999). "Predicting future performance of a shallow steam-zone in the Svartsengi geothermal field, Iceland". *Proceedings of the 24th Workshop on Geothermal Reservoir Engineering*. Stanford University.

Björnsson, G., A. Hjartarson, G. S. Bodvarsson, and B. Stein-grimsson (2003). "Development of a 3-D geothermal reservoir model for the greater Hengill volcano in SW-Iceland". *Proceedings of the TOUGH Symposium 2003*. Lawrence Berkeley National Lab Berkeley, CA.

Borsukiewicz-Gozdur, A. (2010). "Dual-fluid-hybrid power plant co-powered by low-temperature geothermal water". *Geothermics* 39(2), pages 170–176. DOI: 10.1016/j.geothermics.2009.10.004.

Buck, W. R., P. Einarsson, and B. Brandsdóttir (2006). "Tectonic stress and magma chamber size as controls on dike propagation: Constraints from the 1975-1984 Krafla rifting episode". *Journal of Geophysical Research: Solid Earth* 111(B12), n/a–n/a. ISSN: 0148-0227. DOI: 10.1029/2005jb003879.

Buckley, S. J., K. Ringdal, N. Naumann, B. Dolva, T. H. Kurz, J. A. Howell, and T. J. Dewez (2019). "LIME: Software for 3-D visualization, interpretation, and communication of virtual geoscience models". *Geosphere* 15(1), pages 222–235. ISSN: 1553-040X. DOI: 10.1130/ges02002.1.

Burchardt, S. (2008). "New insights into the mechanics of sill emplacement provided by field observations of the Njardvik Sill, Northeast Iceland". *Journal of Volcanology and Geothermal Research* 173(3-4), pages 280–288. ISSN: 0377-0273. DOI: 10.1016/j.jvolgeores.2008.02.009.

Burchardt, S., B. Óskarsson, L. Gustafsson, S. E. Berg, and M. S. Riishuus (2022). "Geology of a Neogene caldera cluster in the Borgarfjörður eystri–Loðmundarfjörður area, Eastern Iceland". *Volcanica* 5(1), pages 133–161. DOI: 10.30909/vol.05.01.133161.

Burchardt, S. and T. R. Walter (2010). "Propagation, linkage, and interaction of caldera ring-faults: comparison between analogue experiments and caldera collapse at Miyakejima, Japan, in 2000". *Bulletin of Volcanology* 72(3), pages 297–308. ISSN: 1432-0819. DOI: 10.1007/s00445-009-0321-7.

Burchardt, S., T. R. Walter, and H. Tuffen (2018). "Growth of a Volcanic Edifice Through Plumbing System Processes—Volcanic Rift Zones, Magmatic Sheet-Intrusion Swarms and Long-Lived Conduits". *Volcanic and Igneous Plumbing Systems*, pages 89–112. DOI: 10.1016/b978-0-12-809749-6.00004-2.

Chanceaux, L. and T. Menand (2014). "Solidification effects on sill formation: An experimental approach". *Earth and Planetary Science Letters* 403, pages 79–88. ISSN: 0012-821X. DOI: 10.1016/j.epsl.2014.06.018.

Cubuk-Sabuncu, Y., K. Jónsdóttir, C. Caudron, T. Lecocq, M. M. Parks, H. Geirsson, and A. Mordret (2021). "Temporal Seismic Velocity Changes During the 2020 Rapid Inflation at Mt. Þorbjörn-Svartsengi, Iceland, Using Seismic Ambient Noise". *Geophysical Research Letters* 48(11). DOI: 10.1029/2020gl092265.

Di Muro, A., F. Schwarzlmuller, U. Kueppers, M. Heap, and D. Dingwell (2021). "Petrophysical characterisation of volcanic ejecta to constrain subsurface lithological heterogeneities: implications for edifice stability at basaltic volcanoes". *Volcanica* 4(1), pages 41–66. ISSN: 2610-3540. DOI: 10.30909/vol.04.01.4166.

Ebigbo, A., P. S. Lang, A. Paluszny, and R. W. Zimmerman (2016). "Inclusion-Based Effective Medium Models for the Permeability of a 3D Fractured Rock Mass". *Transport in Porous Media* 113(1), pages 137–158. ISSN: 1573-1634. DOI: 10.1007/s11242-016-0685-z.

Eggertsson, G., Y. Lavallée, J. Kendrick, and S. Markússon (2020). "Improving fluid flow in geothermal reservoirs by thermal and mechanical stimulation: The case of Krafla volcano, Iceland". *Journal of Volcanology and Geothermal Research* 391, page 106351. ISSN: 0377-0273. DOI: 10.1016/j.jvolgeores.2018.04.008.

Eichelberger, J. (2019). "Planning an International Magma Observatory". *Eos* 100. DOI: 10.1029/2019eo125255.

- (2020). "Distribution and Transport of Thermal Energy within Magma-Hydrothermal Systems". *Geosciences* 10(6), page 212. DOI: 10.3390/geosciences10060212.

Einarsson, P. (1991). "Earthquakes and present-day tectonism in Iceland". *Tectonophysics* 189(1-4), pages 261–279. DOI: 10.1016/0040-1951(91)90501-i.

Elders, W. A., G. Ó. Friðleifsson, R. A. Zierenberg, E. C. Pope, A. K. Mortensen, Á. Guðmundsson, J. B. Lowenstern, N. E. Marks, L. Owens, D. K. Bird, M. Reed, N. J. Olsen, and P. Schiffman (2011). "Origin of a rhyolite that intruded a geothermal well while drilling at the Krafla volcano, Iceland". *Geology* 39(3), pages 231–234. ISSN: 0091-7613. DOI: 10.1130/g31393.1.

<figure><img src="image_30.png" /></figure>
Elders, W. A., G. Ó. Friðleifsson, and B. Pálsson (2014). "Iceland Deep Drilling Project: The first well, IDDP-1, drilled into magma". *Geothermics* 49, page 1. ISSN: 0375-6505. DOI: 10.1016/j.geothermics.2013.08.012.

Eldursi, K., Y. Branquet, L. Guillou-Frottier, and E. Marcoux (2009). "Numerical investigation of transient hydrothermal processes around intrusions: Heat-transfer and fluid-circulation controlled mineralization patterns". *Earth and Planetary Science Letters* 288(1-2), pages 70-83. ISSN: 0012-821X. DOI: 10.1016/j.epsl.2009.09.009.

Farquharson, J., M. J. Heap, N. R. Varley, P. Baud, and T. Reuschlé (2015). "Permeability and porosity relationships of edifice-forming andesites: A combined field and laboratory study". *Journal of Volcanology and Geothermal Research* 297, pages 52-68. ISSN: 0377-0273. DOI: 10.1016/j.jvolgeores.2015.03.016.

Farquharson, J. I., B. Wild, A. R. L. Kushnir, M. J. Heap, P. Baud, and B. Kennedy (2019). "Acid-Induced Dissolution of Andesite: Evolution of Permeability and Strength". *Journal of Geophysical Research: Solid Earth* 124(1), pages 257-273. DOI: 10.1029/2018jb016130.

Foumier, R. O. (1991). "The transition from hydrostatic to greater than hydrostatic fluid pressure in presently active continental hydrothermal systems in crystalline rock". *Geophysical Research Letters* 18(5), pages 955-958. ISSN: 0094-8276. DOI: 10.1029/91gl00966.

Fridleifsson, G. O. and W. A. Elders (2005). "The Iceland Deep Drilling Project: a search for deep unconventional geothermal resources". *Geothermics* 34(3), pages 269-285. ISSN: 0375-6505. DOI: 10.1016/j.geothermics.2004.11.004.

Garden, T. O., I. Chambefort, D. M. Gravley, C. D. Deering, and B. M. Kennedy (2020). "Reconstruction of the fossil hydrothermal system at Lake City caldera, Colorado, U.S.A.: Constraints for caldera-hosted geothermal systems". *Journal of Volcanology and Geothermal Research* 393, page 106794. DOI: 10.1016/j.jvolgeores.2020.106794.

Garden, T. O., D. M. Gravley, B. M. Kennedy, C. Deering, and I. Chambefort (2017). "Controls on hydrothermal fluid flow in caldera-hosted settings: Evidence from Lake City caldera, USA". *Geosphere* 13(6), pages 1993-2016. ISSN: 1553-040X. DOI: 10.1130/ges01506.1.

Gerdes, M. L., L. P. Baumgartner, and M. Person (1998). "Convective fluid flow through heterogeneous country rocks during contact metamorphism". *Journal of Geophysical Research: Solid Earth* 103(B10), pages 23983-24003. ISSN: 0148-0227. DOI: 10.1029/98jb02049.

Ghiorso, M. S. and R. O. Sack (1995). "Chemical mass transfer in magmatic processes IV. A revised and internally consistent thermodynamic model for the interpolation and extrapolation of liquid-solid equilibria in magmatic systems at elevated temperatures and pressures". *Contributions to Mineralogy and Petrology* 119(2-3), pages 197-212. DOI: 10.1007/bf00307281.

Girona, T., V. Realmuto, and P. Lundgren (2021). "Large-scale thermal unrest of volcanoes for years prior to eruption". *Nature Geoscience* 14(4), pages 238-241. DOI: 10.1038/s41561-021-00705-4.

Gunnarsson, G. and E. S. P. Aradóttir (2015). "The Deep Roots of Geothermal Systems in Volcanic Areas: Boundary Conditions and Heat Sources in Reservoir Modeling". *Transport in Porous Media* 108(1), pages 43-59. ISSN: 1573-1634. DOI: 10.1007/s11242-014-0328-1.

Gunnarsson, G., A. Arnaldsson, and A. L. Oddsdóttir (2011). "Model Simulations of the Hengill Area, Southwestern Iceland". *Transport in Porous Media* 90(1), pages 3-22. ISSN: 1573-1634. DOI: 10.1007/s11242-010-9629-1.

Hanson, R. B. and M. D. Barton (1989). "Thermal development of low-pressure metamorphic belts: Results from two-dimensional numerical models". *Journal of Geophysical Research: Solid Earth* 94(B8), pages 10363-10377. ISSN: 0148-0227. DOI: 10.1029/jb094ib08p10363.

Hayba, D. O. and S. E. Ingebritsen (1997). "Multiphase groundwater flow near cooling plutons". *Journal of Geophysical Research: Solid Earth* 102(B6), pages 12235-12252. DOI: 10.1029/97jb00552.

Heap, M. J., D. M. Gravley, B. M. Kennedy, H. A. Gilg, E. Bertolett, and S. L. Barker (2020a). "Quantifying the role of hydrothermal alteration in creating geothermal and epithermal mineral resources: The Ohakuri ignimbrite (Taupo Volcanic Zone, New Zealand)". *Journal of Volcanology and Geothermal Research* 390, page 106703. ISSN: 0377-0273. DOI: 10.1016/j.jvolgeores.2019.106703.

Heap, M. J. and B. M. Kennedy (2016). "Exploring the scale-dependent permeability of fractured andesite". *Earth and Planetary Science Letters* 447, pages 139-150. ISSN: 0012-821X. DOI: 10.1016/j.epsl.2016.05.004.

Heap, M. J., B. M. Kennedy, J. I. Farquharson, J. Ashworth, K. Mayer, M. Letham-Brake, T. Reuschlé, H. A. Gilg, B. Scheu, Y. Lavallée, P. Siratovich, J. Cole, A. D. Jolly, P. Baud, and D. B. Dingwell (2017). "A multidisciplinary approach to quantify the permeability of the Whakaari/White Island volcanic hydrothermal system (Taupo Volcanic Zone, New Zealand)". *Journal of Volcanology and Geothermal Research* 332, pages 88-108. ISSN: 0377-0273. DOI: 10.1016/j.jvolgeores.2016.12.004.

Heap, M. J., A. R. Kushnir, J. Vasseur, F. B. Wadsworth, P. Harlé, P. Baud, B. M. Kennedy, V. R. Troll, and F. M. Deegan (2020b). "The thermal properties of porous andesite". *Journal of Volcanology and Geothermal Research* 398, page 106901. ISSN: 0377-0273. DOI: 10.1016/j.jvolgeores.2020.106901.

Heap, M. J., V. R. Troll, A. R. L. Kushnir, H. A. Gilg, A. S. D. Collinson, F. M. Deegan, H. Darmawan, N. Seraphine, J. Neuberg, and T. R. Walter (2019). "Hydrothermal alteration of andesitic lava domes can lead to explosive volcanic behaviour". *Nature Communications* 10(1). DOI: 10.1038/s41467-019-13102-8.

Heap, M. J., M. Villeneuve, F. Albino, J. I. Farquharson, E. Brothelande, F. Amelung, J.-L. Got, and P. Baud (2020c). "Towards more realistic values of elastic moduli for volcano modelling". *Journal of Volcanology and Geothermal Research* 390, page 106684. ISSN: 0377-0273. DOI: 10.1016/j.jvolgeores.2019.106684.

Heap, M. J. and M. E. Violay (2021). "The mechanical behaviour and failure modes of volcanic rocks: a review". *Bulletin of Volcanology* 83(5). ISSN: 1432-0819. DOI: 10.1007/s00445-021-01447-2.

<figure><img src="image_31.png" /></figure>
Hjartardóttir, Á. R., P. Einarsson, E. Bramham, and T. J. Wright (2012). "The Krafla fissure swarm, Iceland, and its formation by rifting events". *Bulletin of Volcanology* 74(9), pages 2139-2153. ISSN: 1432-0819. DOI: 10.1007/s00445-012-0659-0.

Hobé, A., H. Geirsson, A. Tryggvason, O. Gudmundsson, et al. (2021). "Seismicity Patterns Due to Magma Intrusions Underneath Geothermal Power Plants". *Proceedings World Geothermal Congress 2020+1*. Reykjavik, Iceland, page 1.

Hobé, A., D. Vogler, M. P. Seybold, A. Ebigbo, R. R. Settgast, and M. O. Saar (2018). "Estimating fluid flow rates through fracture networks using combinatorial optimization". *Advances in Water Resources* 122, pages 85-97. ISSN: 0309-1708. DOI: 10.1016/j.advwatres.2018.10.002.

Jónasson, K. (1994). "Rhyolite volcanism in the Krafla central volcano, north-east Iceland". *Bulletin of Volcanology* 56(6-7), pages 516-528. ISSN: 1432-0819. DOI: 10.1007/bf00302832.

Karaoğlu, Ö., M. Bazargan, A. Baba, and J. Browning (2019). "Thermal fluid circulation around the Karliova triple junction: Geochemical features and volcano-tectonic implications (Eastern Turkey)". *Geothermics* 81, pages 168-184. DOI: 10.1016/j.geothermics.2019.05.003.

Karaoğlu, Ö., J. Browning, M. Bazargan, and A. Gudmundsson (2016). "Numerical modelling of triple-junction tectonics at Karliova, Eastern Turkey, with implications for regional magma transport". *Earth and Planetary Science Letters* 452, pages 157-170. DOI: 10.1016/j.epsl.2016.07.037.

Kennedy, B., E. Holohan, J. Stix, D. Gravley, J. Davidson, and J. Cole (2018). "Magma plumbing beneath collapse caldera volcanic systems". *Earth-Science Reviews* 177, pages 404-424. ISSN: 0012-8252. DOI: 10.1016/j.earscirev.2017.12.002.

Kennedy, B. M., M. J. Heap, S. Burchardt, M. Villeneuve, H. Tuffen, H. A. Gilg, J. Davidson, N. Duncan, E. Saubin, E. B. Gestsson, M. Anjomrouz, and P. Butler (2022). "Thermal impact of dykes on ignimbrite and implications for fluid flow compartmentalisation in calderas". *Volcanica* 5(1), pages 75-93. ISSN: 2610-3540. DOI: 10.30909/vol.05.01.7593.

Khoei, A. R., M. Hirmand, M. Vahab, and M. Bazargan (2015). "An enriched FEM technique for modeling hydraulically driven cohesive fracture propagation in impermeable media with frictional natural faults: Numerical and experimental investigations". *International Journal for Numerical Methods in Engineering* 104(6), pages 439-468. DOI: 10.1002/nme.4944.

Larsen, G., K. Grönvold, and S. Thorarinsson (1979). "Volcanic eruption through a geothermal borehole at Námafjall, Iceland". *Nature* 278(5706), pages 707-710. DOI: 10.1038/278707a0.

Lehmann, J., N. Arndt, B. Windley, M.-F. Zhou, C. Y. Wang, and C. Harris (2007). "Field Relationships and Geochemical Constraints on the Emplacement of the Jinchuan Intrusion and its Ni-Cu-PGE Sulfide Deposit, Gansu, China". *Economic Geology* 102(1), pages 75-94. ISSN: 0361-0128. DOI: 10.2113/gsecongeo.102.1.75.

Liotta, D., A. Brogi, G. Ruggieri, V. Rimondi, M. Zucchi, H. M. Helgadóttir, G. Montegrossi, and G. Ó. Friðleifsson (2020). "Fracture analysis, hydrothermal mineralization and fluid pathways in the Neogene Geitafell central volcano: insights for the Krafla active geothermal system, Iceland". *Journal of Volcanology and Geothermal Research* 391, page 106502. DOI: 10.1016/j.jvolgeores.2018.11.023.

Marsh, B. D. (1989). "On Convective Style and Vigor in Sheet-like Magma Chambers". *Journal of Petrology* 30(3), pages 479-530. DOI: 10.1093/petrology/30.3.479.

Mathieu, L., S. Burchardt, V. R. Troll, M. Krumbholz, and A. Delcamp (2015). "Geological constraints on the dynamic emplacement of cone-sheets – The Ardnamurchan cone-sheet swarm, NW Scotland". *Journal of Structural Geology* 80, pages 133-141. ISSN: 0191-8141. DOI: 10.1016/j.jsg.2015.08.012.

Mikofski, M. (2013). *IAPWS*†*F97v0.1*. https://github.com/mikofski/IAPWS_IF97.

Mordensky, S. P., B. M. Kennedy, M. C. Villeneuve, Y. Lavalée, M. K. Reichow, P. A. Wallace, P. A. Siratovich, and D. M. Gravley (2019). "Increasing the Permeability of Hydrothermally Altered Andesite by Transitory Heating". *Geochemistry, Geophysics, Geosystems* 20(11), pages 5251-5269. DOI: 10.1029/2019gc008409.

Mordensky, S., M. Villeneuve, B. Kennedy, M. Heap, D. Gravley, J. Farquharson, and T. Reuschlé (2018). "Physical and mechanical property relationships of a shallow intrusion and volcanic host rock, Pinnacle Ridge, Mt. Ruapehu, New Zealand". *Journal of Volcanology and Geothermal Research* 359, pages 1-20. ISSN: 0377-0273. DOI: 10.1016/j.jvolgeores.2018.05.020.

Mortensen, A., G. Ásgrímur, B. Steingrímsson, G. Axelsson, H. Ármannsson, H. Björnsson, K. Ágústsson, K. Saemundsson, M. Ólafsson, R. Karlsdóttir, S. Halldórsdottir, and T. Hauksson (2015). *The Krafla Geothermal System. Research summary and conceptual model revision*. Technical report. Landsvirkjun LV-2015-098.

Mortensen, A., Þ. Egilson, B. Gautason, S. Árnadóttir, and Á. Guðmundsson (2014). "Stratigraphy, alteration mineralogy, permeability and temperature conditions of well IDDP-1, Krafla, NE-Iceland". *Geothermics* 49, pages 31-41. ISSN: 0375-6505. DOI: 10.1016/j.geothermics.2013.09.013.

Nichols, A. R. L., M. R. Carroll, and Á. Höskuldsson (2002). "Is the Iceland hot spot also wet? Evidence from the water contents of undegassed submarine and subglacial pillow basalts". *Earth and Planetary Science Letters* 202(1), pages 77-87. DOI: 10.1016/s0012-821x(02)00758-6.

Orkustofnun - National Energy Authority of Iceland (2021). *OS-2021-T008-01: Primary Energy Use in Iceland* 1940-2020. https://orkustofnun.is/gogn/Talnaefni/OS-2021-T008-01.pdf. [Data file].

Óskarsson, B. V. (2015). "Volcanological studies of Neogene flood basalt groups in eastern Iceland". PhD thesis. University of Iceland.

Óskarsson, B. V. and M. S. Riishuus (2013). "The mode of emplacement of Neogene flood basalts in Eastern Iceland: Facies architecture and structure of the Hólmar and Grjótá olivine basalt groups". *Journal of Volcanology and Geothermal Research* 267, pages 92-118. ISSN: 0377-0273. DOI: 10.1016/j.jvolgeores.2013.09.010.

<figure><img src="image_32.png" /></figure>
Paquet, F., O. Dauteuil, E. Hallot, and F. Moreau (2007). "Tectonics and magma dynamics coupling in a dyke swarm of Iceland". *Journal of Structural Geology* 29(9), pages 1477–1493. ISSN: 0191-8141. DOI: 10.1016/j.jsg.2007.06.001.

Parry, W. T., J. C. Bellows, J. S. Gallagher, and A. H. Harvey (2000). *ASME International Steam Tables for industrial use: Based on the IAPWS International Formulation 1997 for the thermodynamic Properties of Water and Steam (IAPWS-IF97)*. Volume CRTD-58. American Society of Mechanical Engineers, New York, N.Y. ISBN: 0791801543.

Putirka, K. D. (2008). "Thermometers and Barometers for Volcanic Systems". *Reviews in Mineralogy and Geochemistry* 69(1), pages 61–120. ISSN: 1529-6466. DOI: 10.2138/rmg.2008.69.3.

Reinsch, T., P. Dobson, H. Asanuma, E. Huenges, F. Poletto, and B. Sanjuan (2017). "Utilizing supercritical geothermal systems: a review of past ventures and ongoing research activities". *Geothermal Energy* 5(1). ISSN: 2195-9706. DOI: 10.1186/s40517-017-0075-y.

Robertson, E. C. and D. L. Peck (1974). "Thermal conductivity of vesicular basalt from Hawaii". *Journal of Geophysical Research* 79(32), pages 4875–4888. DOI: 10.1029/jb079i032p04875.

Rooyakkers, S. M., J. Stix, K. Berlo, M. Petrelli, and F. Sigmundsson (2021). "Eruption risks from covert silicic magma bodies". *Geology*. DOI: 10.1130/g48697.1.

Rubio-Maya, C., V. M. Ambríz Díaz, E. Pastor Martínez, and J. M. Belman-Flores (2015). "Cascade utilization of low and medium enthalpy geothermal resources - A review". *Renewable and Sustainable Energy Reviews* 52, pages 689–716. DOI: 10.1016/j.rser.2015.07.162.

Ruz, J., J. Browning, J. Cembrano, P. Iturrieta, M. Gerbault, and G. Sielfeld (2020). "Field observations and numerical models of a Pleistocene-Holocene feeder dyke swarm associated with a fissure complex to the east of the Tatara-San Pedro-Pellado complex, Southern Volcanic Zone, Chile". *Journal of Volcanology and Geothermal Research* 404, page 107033. ISSN: 0377-0273. DOI: 10.1016/j.jvolgeores.2020.107033.

Satou, C., A. Gudmundsson, R. Gertisser, C. B. Ramsey, M. Bazargan, D. M. Pyle, S. Wulf, A. J. Miles, and M. Hardiman (2021). "Eruptive activity of the Santorini Volcano controlled by sea-level rise and fall". *Nature Geoscience* 14(8), pages 586–592. DOI: 10.1038/s41561-021-00783-4.

Saubin, E., B. Kennedy, H. Tuffen, A. R. L. Nichols, M. Villeneuve, I. Bindeman, A. Mortensen, C. Schipper, F. Wadsworth, T. Watson, and R. Zierenberg (2021). "Textural and geochemical window into the IDDP-1 rhyolitic melt, Krafla, Iceland, and its reaction to drilling". *GSA Bulletin* 133(9–10), pages 1815–1830. ISSN: 1943-2674. DOI: 10.1130/b35598.1.

Saubin, E., B. Kennedy, H. Tuffen, M. Villeneuve, J. Davidson, and S. Burchardt (2019). "Comparative field study of shallow rhyolite intrusions in Iceland: Emplacement mechanisms and impact on country rocks". *Journal of Volcanology and Geothermal Research* 388, page 106691. ISSN: 0377-0273. DOI: 10.1016/j.jvolgeores.2019.106691.

Schaefer, L. N., J. E. Kendrick, T. Oommen, Y. Lavallã©e, and G. Chigna (2015). "Geomechanical rock properties of a basaltic volcano". *Frontiers in Earth Science* 3. ISSN: 2296-6463. DOI: 10.3389/feart.2015.00029.

Schopka, H. H., M. T. Gudmundsson, and H. Tuffen (2006). "The formation of Helgafell, southwest Iceland, a monogenetic subglacial hyaloclastite ridge: Sedimentology, hydrology and volcano-ice interaction". *Journal of Volcanology and Geothermal Research* 152(3–4), pages 359–377. ISSN: 0377-0273. DOI: 10.1016/j.jvolgeores.2005.11.010.

Scott, S., T. Driesner, and P. Weis (2015). "Geologic controls on supercritical geothermal resources above magmatic intrusions". *Nature Communications* 6(1). ISSN: 2041-1723. DOI: 10.1038/ncomms8837.

Sigmundsson, F., A. Hooper, S. Hreinsdóttir, K. S. Vogfjörð, B. G. Ófeigsson, E. R. Heimisson, S. Dumont, M. Parks, K. Spaans, G. B. Gudmundsson, V. Drouin, T. Árnadóttir, K. Jónsdóttir, M. T. Gudmundsson, T. Högnadóttir, H. M. Fridriksdóttir, M. Hensch, P. Einarsson, E. Magnússon, S. Samsonov, B. Brandsdóttir, R. S. White, T. Ágústsdóttir, T. Greenfield, R. G. Green, Á. R. Hjartardóttir, R. Pedersen, R. A. Bennett, H. Geirsson, P. C. La Femina, H. Björnsson, F. Pálsson, E. Sturkell, C. J. Bean, M. Möllhoff, A. K. Braiden, and E. P. S. Eibl (2014). "Segmented lateral dyke growth in a rifting event at Bárðarbunga volcanic system, Iceland". *Nature* 517(7533), pages 191–195. ISSN: 1476-4687. DOI: 10.1038/nature14111.

Siratovich, P. A., M. C. Villeneuve, J. W. Cole, B. M. Kennedy, and F. Bégué (2015). "Saturated heating and quenching of three crustal rocks and implications for thermal stimulation of permeability in geothermal reservoirs". *International Journal of Rock Mechanics and Mining Sciences* 80, pages 265–280. ISSN: 1365-1609. DOI: 10.1016/j.ijrmms.2015.09.023.

Stix, J., B. Kennedy, M. Hannington, H. Gibson, R. Fiske, W. Mueller, and J. Franklin (2003). "Caldera-forming processes and the origin of submarine volcanogenic massive sulfide deposits". *Geology* 31(4), page 375. DOI: 10.1130/0091-7613(2003)031<0375:cfpato>2.0.co;2.

Sturkell, E., P. Einarsson, M. J. Roberts, H. Geirsson, M. T. Gudmundsson, F. Sigmundsson, V. Pinel, G. B. Guðmundsson, H. Ólafsson, and R. Stefánsson (2008). "Seismic and geodetic insights into magma accumulation at Katla subglacial volcano, Iceland: 1999 to 2005". *Journal of Geophysical Research* 113(B3). ISSN: 0148-0227. DOI: 10.1029/2006jb004851.

Thordarson, T. and G. Larsen (2007). "Volcanism in Iceland in historical time: Volcano types, eruption styles and eruptive history". *Journal of Geodynamics* 43(1), pages 118–152. ISSN: 0264-3707. DOI: 10.1016/j.jog.2006.09.005.

Vogler, E. D. (2014). "Berggrunnskortlagning Breiðuvíkur á Austfjörððum". PhD thesis. University of Iceland.

Walker, G. P. (1999). "Volcanic rift zones and their intrusion swarms". *Journal of Volcanology and Geothermal Research* 94(1-4), pages 21–34. DOI: 10.1016/s0377-0273(99)00096-7.

<figure><img src="image_33.png" /></figure>
Weaver, J., G. H. Eggertsson, J. E. P. Utley, P. A. Wallace, A. Lamur, J. E. Kendrick, H. Tuffen, S. H. Markússon, Y. Lavallée, and E. Gomez-Rivas (2020). "Thermal Liability of Hyaloclastite in the Krafla Geothermal Reservoir, Iceland: The Impact of Phyllosilicates on Permeability and Rock Strength". *Geofluids* 2020, pages 1–20. ISSN: 1468-8115. DOI: 10.1155/2020/9057193.

Weisenberger, T. B., G. Axelsson, A. Arnaldsson, A. Blischke, F. Óskarsson, H. Ármannsson, H. Blanck, H. M. Helgadóttir, J.-C. C. Berthet, K. Ágústsson, S. S. Gylfadóttir, and V. Guðmundsdóttir (2015). *Revision of the Conceptual Model of the Krafla Geothermal System*. Technical report. Landsvirkjun LV-2015-040.

Wohletz, K. and G. Heiken (1992). *Volcanology and geothermal energy*. Volume 432. University of California Press Berkeley.

Wyering, L., M. Villeneuve, I. Wallis, P. Siratovich, B. Kennedy, D. Gravley, and J. Cant (2014). "Mechanical and physical properties of hydrothermally altered rocks, Taupo Volcanic Zone, New Zealand". *Journal of Volcanology and Geothermal Research* 288, pages 76–93. DOI: 10.1016/j.jvolgeores.2014.10.008.

<figure><img src="image_34.png" /></figure>
APPENDIX A

<table><caption>Table A1: Whole rock composition from Öskarsson [2015] used for the equilibrium test. Water content is based on Nichols et al. [2002].</caption><thead><tr><th>Whole rock sample</th><th>SiO<sub>2</sub></th><th>TiO<sub>2</sub></th><th>Al<sub>2</sub>O<sub>3</sub></th><th>FeO</th><th>MnO</th><th>MgO</th><th>CaO</th><th>Na<sub>2</sub>O</th><th>K<sub>2</sub>O</th><th>P<sub>2</sub>O<sub>5</sub></th><th>H<sub>2</sub>O</th><th>Anhydrous total</th></tr></thead><tbody><tr><td>MSR11-010</td><td>49.5</td><td>2.7</td><td>14.6</td><td>13.8</td><td>0.2</td><td>4.6</td><td>8.5</td><td>3.3</td><td>1.1</td><td>0.5</td><td>0.4</td><td>98.8</td></tr></tbody></table>

APPENDIX B

The numerical model domain is subdivided into 11,000 triangular elements (see examples in Figures B1 and B2). We defined three domains with different mesh resolutions: The finest mesh size ensures the highest resolution in, and immediately around, the intrusion(s) to better calculate the heat equation, the second finest mesh is defined in the domain around and above the intrusions to sufficiently resolve the temperature field and fluid transport caused by the intrusion. The mesh concentration inside the intrusion depends on the minimum element size and maximum element size which is different in single thick intrusions and the thin intrusions.

<figure><img src="image_35.png" /><figcaption>Figure B1: Mesh geometry in the models with a single thick sill.</figcaption></figure>

<figure><img src="image_36.png" /></figure>
<figure><img src="image_37.png" /><figcaption>Figure B2: Mesh geometry in the models with eight thin sills.</figcaption></figure>

APPENDIX C

<figure><img src="image_38.png" /><figcaption>Figure C1: Cross section through Leirfjall mountain based on a 3D reconstruction of the subsurface geology. The proportions of intrusions was calculated in ImageJ.</figcaption></figure>

<figure><img src="image_39.png" /></figure>
## APPENDIX D

The entire model domain is governed by Darcy flow through pores in a non-deforming rock matrix.

The Darcy flow law is given by

$$ \mathbf{u} = -\frac{k}{\mu}(\nabla p + \rho_f g) \quad (D1) $$

where $\mathbf{u}$ is the Darcy velocity vector, $k$ is the permeability, $\mu$ is the viscosity of the pore fluid, $p$ is the pore fluid pressure, $\rho_f$ is the pore fluid density, and $g$ is the gravity vector. For a complete description of the mechanical aspects of fluid flow, the Darcy velocity is coupled with a statement of conservation of mass,

$$ \Phi \frac{\partial \rho_f}{\partial t} + \nabla \cdot (\rho_f \mathbf{u}) = 0 \quad (D2) $$

where $\Phi$ is the porosity, which is generally variable in space but constant in time.

In our model, the fluid flow is entirely driven by density differences due to differences in temperature

and pressure. The evolution of the temperature field is dictated by energy conservation,

$$ (\rho c_p)_{\text{eff}} \frac{\partial T}{\partial t} + \rho_f c_p \mathbf{u} \cdot \nabla T = \nabla \cdot (k'_{\text{eff}} \nabla T) \quad (D3) $$

where $c_p$ is the specific heat, $T$ is temperature, and $k'$ is thermal conductivity. The subscript “eff” refers to effective properties of the combined fluid-solid medium. They are defined as

$$ (\rho c_p)_{\text{eff}} = \Phi \rho_f c_{p,f} + (1 - \Phi) \rho_s c_{p,s} \quad (D4) $$

and

$$ k'_{\text{eff}} = \Phi k'_{f} + (1 - \Phi) k'_{s} \quad (D5) $$

where subscript $f$ stands for fluid properties, and subscript $s$ stands for properties pertaining to the solid rock matrix. Use of the energy equation as stated above implies that the temperature of the pore fluid is everywhere in equilibrium with the surrounding rock matrix. Given the low fluid interstitial velocities and the small spatial scale of the pore spaces, thermal diffusion between fluid and solid is expected to be fast relative to the timescales of the large-scale convective flow, and hence the assumption of local thermal equilibrium is justified even in the regions of relatively fast Darcy flow.

The basaltic dykes and sills are intruded in mostly liquid form. The initial intrusion temperature, which is based on thermobarometry, is $T = 1130$ °C. This temperature is just below the liquidus temperature computed from MELTS [Ghiorso and Sack 1995], which is $T_{\text{liq}} = 1145$ °C. As the enthalpy of fusion is released over the temperature range between the liquidus and the solidus, a common practice is to include the enthalpy of fusion in the heat capacity [e.g. Hanson and Barton 1989]. In most models, we assume a linear solidification rate between liquidus and solidus, which leads to addition of a constant to the specific heat over that temperature range. Following Hanson and Barton [1989], we use for these models $T_{\text{liq}} = 950$ °C and $T_{\text{sol}} = 1150$ °C. The enthalpy of fusion is $L = 375 \text{ kJ kg}^{-1}$ [Hanson and Barton 1989]. The specific heat is thus increased from the standard value of $c_{p0} = 1050 \text{ J kg}^{-1} \text{ K}$ to $c_p = 2925 \text{ J kg}^{-1} \text{ K}$ over the temperature range $950 \text{ °C} < T < 1150 \text{ °C}$. With one model, we test the effect of the assumption of linear solidification by using the melting function $f(T)$ as derived from MELTS. For this model, the specific heat is computed using

$$ c_p = c_{p0} + L \frac{df}{dT} \quad (D6) $$

<figure><img src="image_40.png" /></figure>
<figure><img src="image_41.png" /></figure>

Figure D1: Melting functions versus temperature. Most models use the linear melting function (dashed line). In order to evaluate the sensitivity to this approximation, we also ran one model with the melting function calculated using MELTS (solid line). This melting function was calculated at a pressure of 300 MPa for the composition specified in Table A1. For comparison, we also show the melting function for the same composition, but without water (dash-dotted line); this melting function was not used in the model runs.

<figure><img src="image_42.png" /></figure>

Figure D2: Specific heat as function of temperature. The solid line is based on the melting function $f(T)$ from MELTS using the partly wet composition. The dashed line is based on the linear melting assumption. Note that both models use the same value of the enthalpy of fusion $L$, and hence the areas under the two curves are equal. By comparison with the dry melting curve in Figure D1, it is observed that the sharp peak of the specific heat function at 885 °C as well as the elevated values up to 970 °C are due to the presence of water.

<figure><img src="image_43.png" /></figure>
<figure><img src="image_44.png" /><figcaption>Figure D3: Initial geotherm with $dT/dz = 200$ °C in the upper 1000 m, and leveling out to a constant temperature of $T = 200$ °C below a depth of 1100 m, with a smooth transition in between.</figcaption></figure>

<figure><img src="image_45.png" /><figcaption>Figure D4: Permeability of magmatic intrusions as a function of temperature as used in the numerical models. Above the solidus ($T > 900$ °C), intrusions are assumed to be impermeable, and the permeability is set to a low value of $k = 10^{-22}$ m². Between $T = 500$ °C and $T = 900$ °C, permeability is assumed to decrease linearly with increasing $T$ from a value of $k = 10^{-17}$ m² to a value of $10^{-22}$ m² [Hanson and Barton 1989; Hayba and Ingebritsen 1997]. Between $T = 360$ °C and $T = 500$ °C, permeability is taken to decrease linearly from $k = 10^{-16}$ m² to $k = 10^{-17}$ m² with increasing $T$ due to fracture formation in the intrusions [Foumier 1991; Hayba and Ingebritsen 1997].</figcaption></figure>

<figure><img src="image_46.png" /></figure>
<figure><img src="image_47.png" /></figure>

Figure D5: Density of water as a function of temperature for 3 different linear geotherms. Geotherm 1 starts at the surface (1 atm) and 100 °C and ends at a depth of 1500 m and 600 °C. Geotherm 2 starts at 1 atm and 100 °C and ends at 1500 m depth and 400 °C. Geotherm 3 starts at 1 atm and 0 °C and ends at 1500 m depth and 400 °C. Data from Parry et al. [2000], using the MATLAB functions published by Mikofski [2013].

<figure><img src="image_48.png" /></figure>

Figure D6: Specific heat capacity of water as a function of temperature for the same 3 geotherms as in Figure D5. Data from Parry et al. [2000].

<figure><img src="image_49.png" /></figure>
<figure><img src="image_50.png" /><figcaption>Figure D7: Thermal conductivity of water as a function of temperature for the same 3 geotherms as in Figure D5. Data from Parry et al. [2000].</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>Figure D8: Dynamic viscosity of water as a function of temperature for the same 3 geotherms as in Figure D5. Data from Parry et al. [2000].</figcaption></figure>

<figure><img src="image_52.png" /></figure>
<table><caption>Table D1: Permeability functions derived from empirical values from Eggertsson et al. [2020] assuming hydrostatic conditions. Exponential function equation: $k(z) = a_1 \exp(-a_2z) + a_3$; $z_{Min} = 0$ m, $z_{Max} = 1500$ m.</caption><thead><tr><th></th><th>a<sub>1</sub></th><th>a<sub>2</sub></th><th>a<sub>3</sub></th></tr></thead><tbody><tr><td>High permeability Basalt</td><td>9.6490 × 10<sup>-14</sup></td><td>1.8524 × 10<sup>-4</sup></td><td>-</td></tr><tr><td>Low permeability Basalt</td><td>1.5631 × 10<sup>-15</sup></td><td>1.1302 × 10<sup>-4</sup></td><td>-</td></tr><tr><td>High permeability Ignimbrite</td><td>1.0731 × 10<sup>-13</sup></td><td>1.0806 × 10<sup>-3</sup></td><td>-</td></tr><tr><td>Low permeability Ignimbrite</td><td>2.0234 × 10<sup>-15</sup></td><td>1.1868 × 10<sup>-3</sup></td><td>-</td></tr><tr><td>Fractured fit</td><td>1.1594 × 10<sup>-13</sup></td><td>0.0621</td><td>2.2527 × 10<sup>-13</sup></td></tr><tr><td>Fracture network (super-high permeability)</td><td>8.3731 × 10<sup>-13</sup></td><td>0.0621</td><td>1.6269 × 10<sup>-13</sup></td></tr></tbody></table>

<figure><img src="image_53.png" /><figcaption>Figure D9: Exponential permeability functions shown in Table D1 compared to the experimental permeability as a function of effective pressure data for intact and fractured basalts from Eggertsson et al. [2020]. Conversion to depth was done assuming hydrostatic conditions. The super-high permeability models (fracture network) use a shifted curve corresponding to a fractured basalt (fractured fit). The shift has the curve start at $1 \times 10^{-12}$ at 0 m depth.</figcaption></figure>

<figure><img src="image_54.png" /></figure>