Hyaloclastites commonly form high-quality reservoir rocks in volcanic geothermal provinces. Here, we investigated the effects of confinement due to burial following prolonged accumulation of eruptive products on the physical and mechanical evolution of surficial and subsurface (depths of 70 m, 556 m, and 732 m) hyaloclastites from Krafla volcano, Iceland. Upon loading in a hydrostatic cell, the porosity and permeability of the surficial hyaloclastite decreased linearly with mean effective stress, as pores and cracks closed due to elastic (recoverable) compaction up to 22-24 MPa (equivalent to ~1.3 km depth in the reservoir). Beyond this mean effective stress, denoted as , we observed accelerated porosity and permeability reduction with increasing confinement, as the rock underwent permanent inelastic compaction. In comparison, the porosity and permeability of the subsurface core samples were less sensitive to mean effective stress, decreasing linearly with increasing confinement as the samples compacted elastically within the conditions tested (to 40 MPa). Although the surficial material underwent permanent, destructive compaction, it maintained higher porosity and permeability than the subsurface hyaloclastites throughout the experiments. We constrained the evolution of yield curves of the hyaloclastites, subjected to different effective mean stresses in a triaxial press. Surficial hyaloclastites underwent a brittle-ductile transition at an effective mean stress of ~10.5 MPa, and peak strength (differential stress) reached 13 MPa. When loaded to effective mean stresses of 33 and 40 MPa, the rocks compacted, producing new yield curves with a brittle-ductile transition at ~12.5 and ~19 MPa, respectively, but showed limited strength increase. In comparison, the subsurface samples were found to be much stronger, displaying higher strengths and brittle-ductile transitions at higher effective mean stresses (i.e., 37.5 MPa for 70 m sample, >75 MPa for 556 m, and 68.5 MPa for 732 m) that correspond to their lower porosities and permeabilities. Thus, we conclude that compaction upon burial alone is insufficient to explain the physical and mechanical properties of the subsurface hyaloclastites present in the reservoir at Krafla volcano. Mineralogical alteration, quantified using SEM-EDS, is invoked to explain the further reduction of porosity and increase in strength of the hyaloclastite in the active geothermal system at Krafla.

1. Introduction

Geothermal and hydrothermal systems are typically found in active volcanic environments [14], where fluid convection transfers heat and mass from the high-temperature subsurface (e.g., [5]). As magma underrooted systems can be intermittently volcanically active over long periods of time (i.e., ~1 Ma), it is common for the reservoir rock, hosting high-enthalpy fluids, to be of volcanic origin. The initial geomechanical properties, such as permeability and strength, of these reservoir rocks can be as varied as the style of volcanism from which they are formed [68] and may be susceptible to subsequent changes due to burial, heat flux, and interaction with saturated fluids (e.g., [9]). Thus, the evolution of geothermal systems fed by magmatic bodies is intrinsically linked to the petrological and mechanical evolution of the reservoir rocks.

In mid- to high-latitude provinces, where volcanic activity may commonly be subaqueous or subglacial (e.g., in Iceland, Chile, and New Zealand), substantially increased cooling rates and elevated pressure promote quench-induced fragmentation and suppress exsolution fragmentation [10]. During basaltic eruptions, the products of such activity include highly variable, quench-fragmented glass, termed hyaloclastite (e.g., [1118]). Through time, the glass commonly undergoes extensive alteration, resulting in the generation of a palagonite matrix, dominated by minerals such as smectite and zeolites, which commonly breakdown when subjected to moderately high temperatures of a few hundred degrees (e.g., [1922]). As such, hyaloclastite comprises a time- and temperature-dependent, variably indurated, and heterogeneous assortment of palagonite, hydrated glass [17], lithics, and crystal fragments. The common presence of clay phases in reservoir rocks can be mapped from the surface using electrical resistivity, providing information about the structure of the reservoir [2325].

Hyaloclastites are often highly porous [2628], mechanically weak [2830], and thus highly permeable and susceptible to fluid circulation [27, 31]. As such, they are often targeted as preferred reservoir rocks for freshwater aquifers [32] and for hydrothermal fluid extraction in geothermal energy production [3335]. In active volcanic systems, hyaloclastites are progressively buried by recurrent deposition of eruptive products and intruded by magmatic bodies, such that they are increasingly in contact with, and host to, hydrothermal fluids [25, 26, 31, 36]. Thus, they experience elevated pressures, high temperatures, and corrosive fluids [26, 3740]. Such extreme conditions may promote compaction [27, 4144], precipitation of secondary mineral phases [26, 4547], and variable degrees of alteration [43, 48], modifying the mechanical properties and the permeable-porous network through which fluids circulate. Previous mechanical studies of porous rock compaction [41, 49, 50] have characterised rock strength by evaluating yield curves to identify the stress conditions where permanent inelastic deformation may occur. These investigations have shown that yielding behaviour can vary between different rock types, with many granular materials (e.g., soils and sandstones) typically having elliptical-shaped yield curves [49, 51], whereas volcanic rocks have been reported to exhibit linear yield curves [52, 53]. Despite their importance in geothermal fields, the mechanical properties and yielding behaviour of hyaloclastite remain largely unconstrained. An understanding of how the pore space and mechanical strength evolve during compaction of hyaloclastite are central to our ability to adequately model the evolution of hydrothermal reservoirs for optimising energy production.

Here, we have systematically mapped the physical and mechanical properties of surficial and subsurface hyaloclastites from Krafla volcano, northeast Iceland, which constitute important reservoir rocks in the active hydrothermal system exploited for geothermal energy. The volcano consists of a large caldera that formed at ~100 ka, possibly in two eruptions [54], that has been partly infilled with ignimbrites, lava flows (commonly occurring every 300-1000 years during the Holocene [54]), hyaloclastites, and other fragmental products [25, 54]. Cores and drill cuttings, obtained from extensive geothermal exploration of the Krafla hydrothermal system, have revealed that the upper >1300 m of the reservoir mostly consists of basaltic lavas and hyaloclastites. Below this depth, the reservoir consists of gabbroic intrusions [25, 36], which are locally underrooted by rhyolitic magma at a depth of 2100 m [55]. As observed in other areas of Iceland, the hydrothermal reservoir rocks can be divided into five zones based on temperature-induced alteration: (1) a shallow zone of smectite-zeolite, (2) interlayered smectite-chlorite, (3) chlorite, (4) chlorite-epidote, and (5) epidote-actinolite [25]. Calcite can additionally occur in regions where the rock is at temperatures lower than ~290°C [25]. Alteration and temperature vary considerably across the field and do not increase linearly with depth [24]; for example, previous studies on core samples from borehole KH-6, drilled in 2006, found that the rocks at 556 m depth were relatively unaltered, but the rock samples collected at 732 m depth had experienced a high intensity of alteration [24].

2. Materials and Methods

In this study, we investigate the effect of confinement on the physical and mechanical evolution of surficial and subsurface hyaloclastites from Krafla volcano, northeast Iceland, to assess the process of compaction and the degree to which compaction contributes to permeability evolution during burial in the reservoir. The surficial sample was collected at the southeastern edge of the caldera (65°N 41.067; -16°W 43.089) during a field campaign in August 2015. It is a basaltic hyaloclastite produced during a subglacial eruption, shortly after the formation of the caldera at ~100 ka [25, 54]. Subsurface samples were selected in August 2016 from cores drilled and collected by Landsvirkjun National Power Company of Iceland: hyaloclastites from a depth of 70-76 m were selected from borehole KH-4 (65°N 41.411; 16°W 48.140) drilled in 2006, and hyaloclastites from 556 m and 732 m depth were selected from borehole KH-6 (65°N 42.115; 16°W 48.048) drilled in 2007 [56]. The sample from 556 m was located approximately 1 m from a basaltic dyke; dykes regularly intrude the subsurface hyaloclastites at Krafla. (Note that within a few centimetres from the dyke, the hyaloclastite showed signs of alteration; hence, we selected samples further away.) These core samples are located about 2.5 km from the main region exploited for geothermal energy [25]. During drilling and after completion of the boreholes, temperature was measured at various depth intervals. In the KH-4 well, a temperature of 40°C was recorded at 65 m. In the KH-6 well, temperature measurements were made at different depths and showed slight warming from measurements taken over a 5-day interval. During the last measurements, the temperature was measured to be 50°C at 250 m, increasing rapidly from 75°C to 200°C at 345 m, before decreasing beyond 500 m depth; at 550 m, the temperature was approximately 150°C, and at 735 m depth, 125°C [56]. A sample overview is presented in Table 1.

2.1. Mineralogical and Petrological Analysis

Petrological analysis and phase distribution for samples from the surface, 70 m, 76 m, 556 m, and 732 m were investigated using an optical microscope and an SEM-EDS (QEMSCAN®, Quantitative Evaluation of Minerals by SCANning electron microscopy). Analysis was performed using this automated SEM-EDS system with a 15 kV accelerating voltage and ~5 nA beam current (see [57] for further details). A step size of 20 μm was used to map an area of , per sample to provide an overview of mineral distribution, and a higher resolution 2 μm step size was used to map detailed textures in a area. In each case, two Bruker energy-dispersive X-ray spectrometers (EDS) recorded the discrete secondary X-rays emitted by the sample, which were used to identify and quantify the mineralogy of the thin section by correlation with a reference database comprising known mineral and glass phases. The mineral distribution was summarised numerically by identifying the relative proportions of each mineral in the SEM-EDS images and normalising them to disregard the measured pore space [58]. All samples were thin sectioned perpendicular to the drilling direction.

2.2. Sample Preparation

For experimental purposes, diameter by long cylindrical cores were prepared (~2 : 1 aspect ratio) from the surface material and available well cores from depths of 70 m, 556 m, and 732 m. All core samples were drilled parallel to the drilled well cores. All prepared samples were kept in a drying oven overnight at 75°C and then cooled and stored in a desiccator before any measurements were undertaken.

2.3. Porosity Determination

The porosity of all sample cores was determined using an AccuPyc 1340 Helium Pycnometer from Micromeritics in the Experimental Volcanology and Geothermal Research Laboratory at the University of Liverpool. The device measures the skeletal sample volume (i.e., rock including isolated pores inaccessible to helium) in a 100 cm3 chamber, with an accuracy of ±0.1% of the sample volume. The connected porosity () is then determined via where is the measured skeletal volume and is the volume calculated by the core dimensions.

2.4. Porosity and Permeability Evolution with Pressure

To simulate the impact of hyaloclastite burial, we determined the porosity and permeability changes associated with increasing effective pressure () using a hydrostatic 250 MPa pressure cell from Sanchez Technologies in the Experimental Volcanology and Geothermal Research Laboratory at University of Liverpool. This method was employed for samples with permeability greater than (corresponding to the approximate determination limit of the apparatus). Jacketed, water-saturated samples were loaded in the pressure vessel to the target confining pressure at 5 increments up to 40 MPa; note that 1 km depth would correspond to a confining pressure of approximately 25 MPa assuming a nominal rock density of . During each loading phase, the change in porosity experienced by the compacting sample was determined by measuring the volume of water expelled (±0.05% accuracy) with the sample held at a pore pressure of 1 MPa (see [7, 27]). Subsequently, following >30 minutes of equilibration at the set effective pressure, permeability was measured via the steady-state flow method [59, 60], by exerting a pressure differential of 1 MPa (2 MPa upstream; 1 MPa downstream) and monitoring fluid discharge in the pumps (with ±0.002 ml accuracy). To assess for the need of Klinkenberg [61] or Forchheimer [62] corrections, the pressure gradient was increased and decreased (between 0 and 2 MPa) to ensure the calculated permeability remained constant as flow rate evolved; we found that these corrections were not needed for any of the samples. Following permeability determination, each sample was loaded to the next increment in effective pressure by increasing confinement, whilst monitoring the volume of water expelled from the sample to track pore closure once more and to measure the permeability again.

For samples with a permeability below the detection limit of the apparatus (i.e., ), permeability was quantified using the pulse transient method [63] in a triaxial apparatus in the Rock Deformation Laboratory at the University of Liverpool (see [64]). The sample was fully saturated in water to a pore fluid pressure of 5 MPa. The fluid pressure was then increased by approximately 0.5 MPa on one side of the sample to set a small pressure differential. This pressure differential across the sample then decayed through time, allowing the permeability to be calculated. Once the measurement was completed, the confining pressure was increased to the next increment and the procedure was repeated [63, 64].

2.5. Mechanical Properties

To constrain the elastic limit of the rocks subjected to isotropic loading (), samples were loaded in the hydrostatic cell by incrementally increasing the confining and pore pressures to 46 MPa and 45 MPa, respectively, ensuring that the effective pressure never exceeded 1 MPa. Then, the effective pressure was increased by reducing the pore pressure at a rate of . The volume of water expelled from the sample was monitored in the volumometer of the pump, and the expelled water was used as a proxy for pore volume. This provided the continuous porosity change as a function of effective pressure up to 45 MPa effective pressure, and was defined as the point of negative inflection in porosity-pressure space (after [65]).

The mechanical properties of the samples were further constrained under unconfined and confined conditions in the Experimental Volcanology and Geothermal Research Laboratory at the University of Liverpool. Uniaxial (unconfined) compressive strength (UCS) measurements were conducted using a 5969 Instron uniaxial press (equipped with a 100 kN load cell with a resolution of 100 N and actuator with a testing range of ) where the samples were all brought to failure (defined by a stress drop exceeding 10%) with a strain rate of 10-5 s-1. The measurements were corrected for machine compliance by pressing the pistons directly together under the same loading conditions; this displacement was then subtracted from the displacement measured during the rock tests in real time using the Bluehill® software from Instron. The slope of the linear elastic portion of the stress-strain loading curves was used to calculate Young’s modulus.

Confined conditions were tested using a TRIAX100-300 triaxial press, developed by Sanchez Technologies. The apparatus controls the experimental conditions (up to 300 MPa of axial load and 100 MPa confining pressure) using four Stigma 300 pumps (the pumps operate up to 100 MPa with a resolution of 50 kPa and have a volume of 300 cm3, maximum flow rate of , a resolution of 10-4 cm3, and a volume control and determination accuracy of 0.1%). Here, only two pumps were used as the triaxial tests were done in the absence of pore fluids: the radial confinement was applied using argon in one pump; the axial deformation was controlled using silicon oil in another pump along with a 1.5 kbar Maximator® gas booster (pressure ratio of 1 : 150). The sample assembly consists of the test specimen loaded between alumina cylindrical spacers of 25 mm diameter, jacketed in a 30 cm long Viton® sleeve. Compliance was constrained by axially loading a sample assembly containing a sample of steel (for which the elastic properties had been accurately constrained a priori) to 300 MPa. By subtracting the idealised elastic deformation of the steel during loading, we quantified the compliance as a function of applied axial stress. To test rock samples, a core was placed in the sample assembly and inserted in the press. The confining pressure was slowly increased to a desired testing value using a manual valve; simultaneously, the axial load was automatically maintained 2 MPa higher than the confining pressure to ensure that this was the direction of the greatest principal stress to prevent the axial pistons from receding. During this phase, the samples experienced a small amount of elastic compaction as cracks closed (e.g., [66]).

To determine the strength of the materials, yield curves were plotted in - space, where is the effective mean stress () and is the differential stress . Yield curves were mapped following the procedure of Bedford et al. [41] where a sample is hydrostatically loaded to a given confining pressure (note there was no pore fluid pressure () in these tests), before an axial load was applied () in order to subject the sample to a differential stress. During axial loading, the sample deformed elastically with a quasilinear stress-strain relationship. The stress build-up was monitored until a deviation from linear elastic loading was observed, marking the onset of yield (i.e., permanent inelastic strain), and the axial load was immediately reduced back to 2 MPa above confining pressure to ensure the sample did not accumulate inelastic damage. The and values at the deviation from linear loading were recorded as the yield point, and the same sample was then taken to different confining pressures and the axial loading procedure was repeated at each pressure increment in order to map out the complete yield curve in - space. This procedure is useful when the available material is limited, as is the case when using recovered subsurface core samples, as an entire yield curve can be reconstructed using only one sample. The yield curve intersects the -axis at the hydrostatic yield point (i.e., no differential stress), typically referred to as ; beyond this point, the rock undergoes permanent inelastic compaction and pore collapse [65].

3. Results

3.1. Petrological and Mineralogical Signatures

The surficial hyaloclastite and the subsurface hyaloclastites cored from different depths in the reservoir exhibit contrasting appearance in hand specimen (Figure 1). The surface hyaloclastite is light brown with ~10% dark basaltic scoria fragments, up to 10 mm in size (Figure 1(a)). At 70 m depth, the hyaloclastite is darker compared to the surface samples and consists of ~35% basalt fragments (Figure 1(b)) similar in size and geometry to the surface hyaloclastite. The sample from 556 m depth (Figure 1(c)) is denser than the shallower samples and contains a brown-green matrix with smaller, dark, poorly defined clasts. The sample from 732 m depth has a grey-greenish colour, with larger basalt clasts than the 556 m sample and similar density (Figure 1(d)).

Using optical microscopy, we examine the petrological characteristics of the hyaloclastites further (Figure 2) and supplement this with SEM-EDS (QEMSCAN®) analysis at two step sizes: at 20 µm to explore the large scale distribution of phases across the samples (Figure 3) and at 2 µm to explore phase distribution in the matrix in detail (Figure 4). In addition to the samples from the surface, 70 m, 556 m, and 732 m, an additional sample from 76 m depth (for which we only have a thin section) was examined to explore textural distinctions across short distances. For the surficial sample, the dark basaltic scoria clasts identified in hand specimen are subrounded to subangular vesicular glassy clasts (Figures 2(a)–2(c)). The glassy clasts are within a heterogeneous matrix of glass fragments, crystal fragments, and very fine-grained smectite and zeolite (Figures 2(a)–2(c), 3(a), 3(b), 4(a), and 4(b)). The porosity consists of rounded vesicles within glassy clasts and fragments, as well as more irregular-shaped pores flanking the margins of grains (Figures 2(a), 3(c), and 4(c)). At 70 m depth, the vesicular, glassy clasts are larger and more rounded, and the glass appears more pristine, with the exception of a few oxidised fragments (Figures 2(d)–2(f)). The matrix of the hyaloclastite is darker in plane polarised light (PPL) compared to the surface samples (Figure 2) and lacks zeolite (Figures 3(d), 3(e), 4(d), and 4(e)). Porosity occurs as round vesicles within individual clasts and large, more angular, connected patches between clasts (Figures 3(f) and 4(f)). The sample from 76 m is texturally similar to the surface and 70 m samples, with subrounded glassy clasts of intermediate size (compared to surface and 70 m samples); oxidised fragments are more common and lithic fragments are rare (Figures 2(g)–2(i)). Zeolite is lacking in the matrix, as in the 70 m sample, and calcite is present (Figures 3(g), 3(h), 4(g), and 4(h)). Pore space is more heterogeneously distributed, with smaller maximum size than in the shallower samples, but still consists of round intraclast vesicles and pore space along the margins of fragments (Figures 3(i) and 4(i)). The sample from 556 m depth is red-brown in PPL and has poorly defined clasts, and it is the only sample where veins are observed (Figures 2(j)–2(l)). The glassy clasts and matrix are cut by several mm-long veins that are zeolite rich (Figures 3(j) and 3(k)). The sample is interspersed with smectite and zeolite, which occasionally forms enriched patches, and has a finer matrix compared to the shallower samples (Figures 4(j) and 4(k)). Porosity is limited to fine fractures (Figures 3(l) and 4(l)). The sample from 732 m depth is also red-brown in PPL and has larger, more well-defined basalt clasts than the 556 m sample (Figures 2(m)–2(o)). Zeolite is distributed within the glassy fragments and infills some spherical pores, whilst smectite and actinolite are found throughout the matrix (Figures 3(m), 3(n), 4(m), and 4(n)). Porosity is limited to fine fractures and traces around clasts (Figures 3(o) and 4(o)).

The phase abundance was assessed quantitatively from the 20 µm SEM-EDS data (Table 2, Figure 5). Several differences are noted across the sample suite, yet there is no definitive systematic change in a single phase’s content with depth. Glass is the most abundant phase at all depths (47.5–69.5%), being present as very fine fragments through to clasts of >5 mm, followed by smectite (17.9–37.3%) which dominates the matrix and infills pores in the glassy clasts. The presence of zeolite is noted at the surface, but is absent from 70 and 76 m samples, and again is present at higher abundance at greater depth. Glass and smectite contents increase and then decrease with depth, making way for the zeolite, which is present in the groundmass, as infilling material of round pores in vesicular glass fragments, and in the margins of glassy grains. Calcite, actinolite, and pyrite all make appearances at depth. Calcite forms isolated patches within the matrix from 70 m depth onwards, as does pyrite but at much lower abundance in the 556 m and 732 m samples, whilst actinolite is distributed throughout the matrix in the 556 and 732 m samples.

3.2. Porosity and Permeability Evolution with Pressure

The hyaloclastites are increasingly less porous with burial depth within the geothermal system, with porosity reduced from 39.7% in the surface samples to 22.1% in the 70 m drill core sample and reaching the lowest values of 12.5% and 13.3% in the 556 m and 732 m samples, respectively (Table 3). The porosity range of the 556 m and 732 m samples overlaps, from 11.8 to 13.8% at 556 m and 12.9 to 13.9% at 732 m. As an illustrative example, the permeability at an effective pressure of 4 MPa decreases from at the surface to at 70 m depth as porosity is almost halved and reduces further to a minimum of at 556 m in the highest density and lowest porosity sample, stabilising to at 732 m (Table 3).

To simulate burial conditions, we subjected the shallow (surface and 70 m) hyaloclastite samples to isotropic loading (keeping ) to observe compaction. Samples were initially loaded to high confining pressure (~40 MPa) and pore pressure; then, the pore pressure was gradually reduced to increase the effective pressure whilst continuously monitoring the pore volume. Following an initial consolidation of the sample and the assembly (<2 MPa), the porosity of the surficial sample decreased quasilinearly until an effective pressure of 22-24 MPa was reached; above which, a greater rate of porosity decrease with increasing effective pressure was observed (a steepening of the slope; Figure 6). In contrast, the sample from 70 m depth compacted linearly as effective pressure increased up to 40 MPa.

In a separate run in which the pore volume and permeability of the sample were evaluated at different, noncontinuous pressure increments, we again observed an increase in the reduction rate of both porosity and permeability with effective pressure between the 20.2 MPa and 25.9 MPa measurements for the surface sample. Thus, the elastic limit, , may be constrained at 22-24 MPa for the surficial hyaloclastite. Following the same procedure, the porosity and permeability of the subsurface 70 m hyaloclastite both decreased linearly with hydrostatic pressure, demonstrating that was not reached up to an effective pressure of 40 MPa (Figure 6). The two ways of measuring pore compaction with increasing effective pressure yield similar quantitative changes; compaction is more significant in the more porous surface sample, which reduced by ~8% whilst the lower porosity 70 m samples compacted by ~5% (Figure 6).

The sensitivity of permeability reduction to increasing effective pressure also appears to be controlled by the initial porosity and permeability of the samples. The initially most porous, permeable samples from the surface have the largest decrease in permeability of more than 2 orders of magnitude as effective pressure is increased (Figure 6). The mid-porosity 70 m samples have a slightly smaller permeability reduction across the same range of effective pressures. The much lower porosity samples from 556 and 732 m have low initial permeabilities which are much less sensitive to effective pressure than the shallow samples.

3.3. Mechanical Behaviour of Buried Hyaloclastite

To understand more about their mechanical fingerprint, systematic, repetitive axial loading of two samples of surficial hyaloclastite was conducted (see Supplementary Figure 2), following the procedure of [41], providing reconstructions of the elliptical yield curve. The two yield curves generated for the surficial hyaloclastite were similar in magnitude and exhibited comparable peaks, marking the transition from the brittle regime (where materials strengthen with pressure and fail via localised deformation) to the ductile regime (where materials weaken with pressure and compact via pervasive deformation). The peaks, termed the critical effective mean stress (critical ), occurred at 9.2-12.0 MPa, corresponding to a peak strength (differential stress, ) of ~13 MPa (Figure 7(a)). Where the sample can no longer withhold any shear stress (), the curves (black circles and red triangles; Figure 7(a)) intersect the effective mean stress () axis at a pressure of ~22 MPa, marking .

Two samples of surficial hyaloclastite were then compacted by increasing the effective mean stress beyond , effectively extending to 33 MPa (black boxes; Figure 7(a)) and 40 MPa (blue triangles; Figure 7(a)). The resultant yield curves (also mapped via repetitive loading) achieved similar peak strengths of 14 to 15 MPa but the curves were elongated, with shifts in the brittle-ductile transition (critical ) to a higher effective mean stress of ~12.5 MPa for the sample compacted to 33 MPa and ~19 MPa for the sample compacted to 40 MPa (Figure 7(a)). We further processed the data by normalising each curve against its value, by dividing the effective mean stress () and the differential stress () at each step by the corresponding value for that curve. This was done to allow a direct comparison between different yield curves (after [41, 67]). The normalised data indicate that lowers with the degree of compaction (i.e., with its modified value; Figure 7(b)).

The strength of the subsurface hyaloclastites collected from boreholes is greater than the surface samples, and the yield curves obtained are thus much greater in magnitude (Figure 7(c)). The magnitude of the yield curves increases with decreasing porosity; the sample from 70 m depth exhibited the brittle-ductile transition (critical ) at an effective mean stress of ~37.5 MPa, the lowest porosity sample from 556 m did not cross the brittle-ductile transition within the pressure conditions tested (up to 75 MPa), and the deepest sample from 732 m showed brittle-ductile transition at an effective mean stress of ~68.5 MPa (Figure 7(c)). The pressure for inelastic compaction () of these hyaloclastites was not met during testing up to >43 MPa for the 70 m sample and >75 MPa for the 556 and 732 m samples.

3.4. Physical Controls on Mechanical Characteristics

Porosity of the hyaloclastites decreases with increasing depth within the geothermal system, with those at 556 m and 732 m having narrowly overlapping ranges (Table 3). As porosity is considered the primary control on strength, we compare the mechanical characteristics of the samples from the surface, 70 m, 556 m, and 732 m to porosity (Figure 8). Decrease in porosity corresponds to a decrease in permeability (Figure 8(a)) and an increase in UCS from 5.4 MPa in the surface samples, to 10.3 MPa in the 70 m drill core sample, 37.1 MPa at 556 m, and 40.0 MPa for the deepest samples from 732 m (Figure 8(b)). Similarly, Young’s modulus increases by over an order of magnitude from 0.8 GPa at the surface, to 1.4 GPa at 70 m, 8.6 GPa at 556 m, and 13.1 GPa at 732 m (Figure 8(c)). Although was not achieved in our tests on the buried samples from the drill core, the yield curves of the buried samples have higher magnitude, indicating an increase in from that of the surface hyaloclastite at 22 MPa. The transition between brittle and ductile behaviour, termed the critical effective mean stress, was however observed at the experimental effective pressures tested (Figure 8(d)). With decreasing porosity, the brittle-ductile transition was shifted from an effective mean stress of 10.5 MPa at the surface to 37.5 for the 70 m sample and 68.5 MPa at 732 m, whilst the slightly lower porosity 556 m sample did not meet the critical value at 75 MPa (Table 3).

4. Interpretation and Discussion

By replicating the stress conditions at various depths in the geothermal reservoir, we assess and compare the compaction response of surficial and subsurface hyaloclastites during compression. Compaction may be associated with multiple phenomena (Figure 9): from an increase in effective pressure which can be caused by either local increases in external applied stress (e.g., due to burial from subsequent deposits and lava emplacement or from magma intrusion) or local decreases in pore pressure (e.g., if fluids are drained, excessively extracted, or if the fluid density fluctuates). During simulated burial, the surface hyaloclastite transitions from elastic to inelastic compaction at an effective pressure of 22 MPa (i.e., ; Figure 2), corresponding to a depth of 1.3 km, assuming that the top part of the reservoir is made of layers of basalt and hyaloclastite with a nominal rock density of and a fluid density of (after [68]). Thus, it is likely that some hyaloclastites in the geothermal system at Krafla (logged at depths of at least 1362 m in IDDP-1; [36]) have undergone inelastic compaction at pressures exceeding due to burial. Moreover, vapour-rich hot zones within the reservoir may drop fluid density to below [69], causing locally higher effective pressures that could push the hyaloclastite to . Finally, dykes and sills are abundant throughout the stratigraphy which may have increased local stresses beyond . At higher pressures, i.e., beyond , the rock compacts, which causes the resultant yield curve of the material to widen, pushing the apparent to higher effective pressures without substantially increasing the strength (Figure 7(a)), resulting in a lower ratio (Figure 7(b)).

The samples retrieved from depth display much lower permeability (Figure 8(a)), as anticipated by the reduction in pore volume. The porosity and permeability reduction during loading (increasing effective pressure) is more significant for the more porous, permeable samples from the surface and from 70 m depth (Figure 6), suggesting that initial confinement has the most significant impact on the ability for fluid flow through the materials, as has been noted in other porous, permeable rocks [27, 70]. The samples also show increasing compressive strength and Young’s modulus with burial depth (Table 3, Figures 8(b) and 8(c)) and have a greater strength at a given effective pressure (Figure 7(c)). The subsurface samples have higher magnitude yield curves, with the brittle to ductile transition (critical ) occurring at significantly higher pressures than the surface samples. Using the aforementioned rock density of and fluid density of (after [68]), the overburden would be ~1.2 MPa for the 70 m deep sample, ~9.5 MPa for the 556 m sample, and ~12.5 MPa for the 732 m sample, which places them in the brittle regime at reservoir conditions (Figure 7(c)). Thus, if high differential stresses were to accumulate, they would cause dilatant rupture. The reservoir is located in a divergent, extensional tectonic setting impacted by recurring volcanic activity (e.g., [71, 72]), with a highly varied stress field [73]; as such, brittle failure is not unlikely and would locally enhance fluid circulation (e.g., [7, 27, 74]) within the reservoir.

for the reservoir hyaloclastites would occur at pressures much greater than those likely experienced at the depth at which they were sampled at Krafla. Yet, the reservoir hyaloclastites are expected to have originally exhibited (upon their formation) similar physical and mechanical attributes to those of the surficial hyaloclastites studied herein; so, considering that the surficial hyaloclastite may undergo ductile (compactant) deformation at as little as ~10.5 MPa and that can be exceeded at ~22 MPa, we anticipate that the reservoir hyaloclastites sampled at 556 m and 732 m have likely suffered some degree of inelastic compaction and even exceeded upon burial in the system. Such a mechanism however fails to recreate all the characteristics of the hyaloclastites forming the geothermal reservoir (Figure 7), indicating that the porous permeable network within these rocks has been additionally modified.

Beyond mechanical compaction (Figure 9), several factors may enhance closure of porosity and strengthening with burial, including temperature-induced weakening (e.g., from local magma intrusion) that lowers and enables more complete compaction (see [75], this volume) and interaction with hydrothermal fluids. We observe that the deep hyaloclastites show signs of reactions induced by elevated temperatures (Figures 24, Table 2); 556 m and 732 m samples are altered to blue-green in hand specimen and red-brown in PPL (Figures 1 and 2), and pore infill by secondary mineral precipitation (see Figures 25). In particular, we note that the zeolite fraction initially decreases from 6.9% at the surface to 0.1-0.2% at 70-76 m and then increases to 17.4% and 13.6% in the 556 m and 732 m samples, hosted in pore space, fractures, and in the margins of glass fragments (Figures 3 and 4; Table 2). In addition, the sample from 556 m contains zeolite-rich veins, which, absent in all other samples, may relate to the close proximity (~1 m) of the sample to a dyke. Dykes commonly intrude the hyaloclastites around Krafla. We postulate that zeolite replacement of the glass is significant after burial at low pressure, low temperature. Comparison of the mineralogical assemblage of the hyaloclastites with the alteration chart of Thien et al. [46] helps in the assessment of the conditions of alteration. The surficial zeolite-bearing hyaloclastite indicates a shallow, low-temperature alteration (<50°C and <5 MPa pressure), whereas the 70-76 m depth hyaloclastites contain no zeolites, which suggests that they were altered through interaction with volcanic gases and by progressive interaction with meteoric fluids at low pressure [46]; potentially, dehydration in the near-intrusion 70-76 m samples could explain their absence (cf. [75]).

The increase in smectite with depth, followed by the subsequent reduction (Figure 5), may indicate the base of the low-temperature alteration zone (where smectite is thermodynamically stable) at around 556–732 m. This is often referred to the clay cap within geothermal reservoirs, which is commonly mapped with electrical soundings (e.g., [2325, 76]). Temperature measurements of 120-200°C within the borehole [24, 56] and temperature reversal in borehole KH-6 suggest it is likely that the reduction of smectite in the 536 m and 732 m samples is related to an older temperature profile with a higher geothermal gradient. We also note that actinolite becomes pervasive in the matrix in samples from 556 m and 732 m and that calcite increases from 0% at the surface to 0.9% at 70 m, 5% at 76 m, 1.1% at 556 m, and 2.5% at 732 m. Higher amounts of calcite filling vesicles in the deeper hyaloclastites indicates precipitation from volcanic fluids at higher pressure (25 MPa) and temperature (250°C), which also suggests the rocks may have previously been at a slightly elevated temperature compared to the present [46].

In combination with burial-induced compaction, modification of the pore space by infilling contributes to the low porosity and permeability and increased strength of the subsurface hyaloclastites. Certain aspects of the deformation and alteration history may be gleaned by textural examination. If minerals are restricted to pores, this suggests conditions were such that minerals could precipitate ahead of deformation, which would serve to strengthen the rock mass and ensure it remains coherent; on the other hand, if fractures are present and infilled, it may suggest that stress conditions promoting deformation were reached ahead of those causing reactions. Figure 9 shows how some hyaloclastites may undergo compaction before alteration, whereas others may first undergo alteration, which may (or may not) strengthen the rock sufficiently to prevent compaction, and it shows how in some areas rocks may avoid compaction by protection from adjacent strong lithologies.

Considering that hyaloclastites are variably porous [2628], it is interesting to note that the model of Thien et al. [46] suggests that a reduction in porosity slows the alteration process. This introduces further complexity into the understanding of the hyaloclastite’s evolutionary history, including the respective timing of progressive weakening via leaching by fluids (e.g., [46]) or thermal destabilisation [75] versus strengthening via compaction and secondary mineral precipitation and how these processes may differ as a function of depth. Higher spatial resolution measurements of mineralogy show high variability throughout the borehole [24]. This may in part be due to crosscutting of the hyaloclastite by basaltic dykes (for example, as seen within 1 m of the 556 m sample) which are common and could locally influence alteration products by elevating temperature, and/or it may be due to heterogeneities in the initial hyaloclastite deposits (which, for example, contain differing grain sizes, lithic and crystal contents); thus, interpretations of alteration should be made with caution when considering a limited sample suite. Pressure and temperature conditions are in constant flux in geothermal systems, especially in such active volcanic systems, and even small fluctuations can change the mechanical response of the materials and reactions taking place, highlighting the need for extensive sampling and high-resolution modelling within these systems.

5. Conclusions

This experimental study investigated the mineralogical, physical, and mechanical evolution of hyaloclastite upon burial in the active hydrothermal system at Krafla volcano. During the burial of fresh, surficial hyaloclastite in a reservoir, local pressures will increase causing the physical properties (porosity and permeability) of the material to evolve upon exceeding the elastic limit (), which prompts collapse of the porous network, resulting in reduced porosity and permeability. Pore collapse results in a modification of the yielding behaviour of the rock with a shift in the brittle-ductile transition (critical ) to a higher effective mean stress. After increasing the effective pressure of the surficial hyaloclastite beyond , it is observed that compaction alone does not recreate the physical and mechanical properties of subsurface hyaloclastites, as strength is not correspondingly increased ( is reduced). The subsurface hyaloclastite samples from the reservoir exhibit a progressive enhancement of strength and reduction in porosity and permeability with burial depth. The yield curves of the subsurface samples differ significantly from those produced by compaction of fresh surficial hyaloclastite. Subsurface hyaloclastites do not achieve within the anticipated reservoir conditions at depths of up to 1362 m (logged at IDDP-1; [36]) and remain within the elastic regime. Yet, the buried hyaloclastite samples’ origin as porous surface deposits suggests they may have already experienced inelastic compaction at pressures exceeding . The failure to reproduce the physical and mechanical characteristics of the subsurface hyaloclastites via mechanical compaction of the surface hyaloclastites beyond leads to the conclusion that in isolation, burial-induced compaction is insufficient to generate reservoir-hosted hyaloclastites. We invoke the additional importance of mineralogical alteration and precipitation from hydrothermal fluids occurring at high temperature, which modified the porous permeable network and led to strengthening of the rock, which may have occurred before, during, and after compaction. Mineralogical, physical, and mechanical processes are in constant competition during the evolution of rocks within geothermal systems; small spatial and temporal fluctuations in the local pressure-temperature environment will dictate the ability for fluid to flow and the potential for energy extraction; thus, detailed studies of such processes are required to maximise the energy potential of geothermal systems such as Krafla.

Data Availability

All data is presented in the manuscript and the supplementary material file. Any data not contained herein is available upon request to the authors.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this article.


We acknowledge SEM-EDS technical support of FEI Company, Hillsboro, Oregon, USA. This study has been supported by Landsvirkjun National Power Company of Iceland and financially supported by Landsvirkjun’s Energy Research fund, grants NÝR-17–2015, NÝR-16–2016, and NÝR-20–2017. This research was further supported by a scholarship from the Institute for Risk and Uncertainty at the University of Liverpool, a Starting Grant from the European Research Council (ERC) to Y. Lavallée on Strain Localisation in Magma (SLiM, no. 306488), an Early Career Fellowship of the Leverhulme Trust granted to J.E. Kendrick (ECF-2016-325), and Research Fellowship of the Leverhulme Trust granted to Y. Lavallée (RF-2019-526\4).

Supplementary Materials

The supplementary material contains two supplementary figures of mechanical data. Figure S1: the UCS data showing the stress-strain loading paths for each of the samples tested within the uniaxial press, where the sample is loaded until a stress drop of >10% is observed. Strength increases with increasing sampling depth within the reservoir, and the slope of the stress strain curves is steeper with increasing sampling depth, corresponding to a lower porosity. Figure S2: the raw loading data for each sample, showing the loading paths for each of the samples tested within the triaxial apparatus. The samples are loaded to the target confining pressure and then axially stressed until they exhibit their elastic limit () before removing the axial load, reducing the effective mean stress again. However, if the sample is loaded past (the surface samples), then the effective mean stress is reduced slightly before the axial load is increased on the sample. (Supplementary Materials)