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

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-inﬁll 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, Kraﬂa volcano. These data inform 2D ﬁnite 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


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 hightemperature 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.
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.Mordensky et al. 2018;Heap et al. 2020a;Kennedy et al. 2022].Typical intrusion types in calderas include cryptodomes/laccoliths, 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 Northeast 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].

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.

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).

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.Ármnannsson 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 [Ármnannsson et al. 1987;Eggertsson et al. 2020].Of critical relevance to this project is that, despite the recent basaltic volcanism, there is now ample 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].

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 geothermometry, and model rock properties are estimated from laboratory permeability experiments performed on rocks from Krafla by Eggertsson et al. [2020].

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 Presses universitaires de �rasbourg 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.

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 2 , Al 2 O 3 , 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: 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  [Ab-An] = 0.28±0.11if  > 1050 °C but  [Ab-An] = 0.1±0.05if  < 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 2 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 2 O saturation [cf.Schopka et al. 2006].

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 9 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).

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- ‡ https://www.comsol.com Presses universitaires de �rasbourg 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 agglomeratefilled 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 hostrock 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.

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 hostrocks.The coarser-grained sedimentary rocks contain plagio-  Presses universitaires de �rasbourg clase with both needle-like and tabular shapes in a matrix of indistinguishable material.

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 An 37 to An 74 with an average of An 64 (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 Wo 40 En 39 Fs 21 , with maximum and minimum values of Wo 42 En 47 Fs 29 and Wo 32 En 29 Fs 15 , 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.

Geothermometry
According to the  [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.8kbar).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].

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.

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 mwide 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 mthick 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).

Model temperatures
The initial temperature of the host-rock was set to follow a vertical temperature gradient of 200 • C km −1 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 9 s (31.7 years) before adding the intrusion(s).
The intrusions are introduced instantaneously with an initially uniform temperature of  = 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.  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].
Temperature Since we consider the intrusions to be impermeable at temperatures above the solidus, we set permeability to 1 × 10 −22 m 2 at  > 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].

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 [Siratovich 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) 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 × 10 −14 m 2 , and low-permeability rock with a permeability of 9.8 × 10 −16 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 perme-Presses universitaires de �rasbourg 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 possible that the geothermal flow is strongly controlled by a single, 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 observe any large veins, we deem it most appropriate to model fractured rock units in terms of isotropic Darcy flow with increased permeability.
The combination of the four intrusion geometries and three host-rock permeabilities results in a total of 12 numerical models, as listed in Table 3.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 transport is dominant in all models with low-permeability hostrocks, while heat from the intrusions in the models with high and super-high permeability host-rocks triggers vigorous convection 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 intrusions create a narrow thermal aureole by heating their immediate 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 intrusions, 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).

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

Presses universitaires de �rasbourg
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 −3 .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].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.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].
Presses universitaires de �rasbourg  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).
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-Presses universitaires de �rasbourg 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 hostrocks, 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 threedimensional 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 emplacemet 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 Presses universitaires de �rasbourg 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 [Árnason 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 flowrates 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); 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]

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.5(2): 477-507.https://doi.org/10.30909/vol.05.02.477507Presses universitaires de �rasbourg 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.Presses universitaires de �rasbourg   [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 2 to k = 10 −17 m 2 with increasing T due to fracture formation in the intrusions [Foumier 1991;Hayba and Ingebritsen 1997].

Figure 1 :
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].

Figure 2 :Figure 3 :Figure 4 :
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 5 :
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.

Figure 6 :
Figure 6: Results from equilibrium test and thermobarometry for plagioclase mineral in three representative intrusions (dykes and sills) in Leirfjall.

Figure 7 :
Figure7: Sketches of the model geometries with intrusions in red and host-rock in grey.Yellow circles mark the locations of temperature measurements in Figure10.Note that the drawings are not exactly to scale to allow visibility of small features.1 atm = 101,325 Pa.

Figure 8 :
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 Eggertsson et al. [2020].See also Appendix D, Figure D3, D4, D5, D6, and D7.

Figure 9 :
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 Table3for a list of all models.See text for description.

Figure 10 :
Figure10: Cooling of the intrusions in all models, measured 20 m from the tip of the (top) sill and (middle) dyke (see locations in Figure7).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 Table3for 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.

Figure 11 :
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).
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 2 sequestration.

Figure B2 :
Figure B2: Mesh geometry in the models with eight thin sills.

Figure D1 :
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 TableA1.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 D2 :
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 D3 :
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.

Figure D4 :
FigureD4: 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 2 .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 2 to a value of 10 −22 m 2[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 2 to k = 10 −17 m 2 with increasing T due to fracture formation in the intrusions[Foumier 1991;Hayba and Ingebritsen 1997].

Figure D5 :
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 D6 :
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 D7 :
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].

Figure D8 :
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].

Figure D9 :
Figure D9: Exponential permeability functions shown in TableD1compared to the experimental permeability as a function of effective pressure data for intact and fractured basalts fromEggertsson 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 × 10 −12 at 0 m depth.

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 superhigh permeability).
or induced via mechanical hydrofracturing [e.g.Siratovich et al. 2015] or chemical stimulation [e.g.Farquharson Presses universitaires de �rasbourg et al. 2019] could herald efficient heat extraction from small intrusions in such host-rocks.