Combined thermo-mechanical jointing experiments

To test the thermo-mechanics of columnar joints we developed a novel experimental setup that allows us to directly observe fracturing in cooling lavas (Supplementary Fig. 1). In the experiment, a tensile fracture is induced by cooling a cylindrical sample of rock from its solidus temperature while locking the ends of the cylinder in position, initially imposing 0 or 5 MPa of uniaxial compressive stress, perpendicular to the experimentally generated fracture plane (thus simulating normal stress at depth). This setup is designed to simulate elastic stress accumulation by thermal contraction between the static centres of two colonnades. Experiments were conducted on a typical, micro-crystalline basaltic lava from a jointed body at Seljavellir, at the base of Eyjafjallajökull, Iceland (Fig. 2a). The colonnades exhibit quasi-hexagonal fracture patterns, jointed into columns ranging between 30 and 130 cm across. The approximately regular spacing of striae, which sometimes pinch in and out laterally (Fig. 2a), scales with the column width16,17,24 (Fig. 2b). The basalt consists of plagioclase, olivine, occasional pyroxenes and iron oxide crystals (Supplementary Figs. 2 and 3), set in a micro-crystalline groundmass of the same mineralogy, which hosts no interstitial glass (Supplementary Fig. 3). The basalt has a solidus temperature of 980 °C (at 1 bar and fO 2 between nickel-nickel oxide (NNO) and quartz-fayalite-magnetite (QFM) oxygen buffers) estimated by MELTS30 on the basis of the bulk rock chemistry (Supplementary Table 1); this was supported visually by the onset of melting of iron oxides at the starting temperature of the experiments. In the columnar jointing experiments, cooling resulted in tensile stress build-up from the solidus temperature down to temperatures of 890–840 °C, regardless of the imposed cooling rate and initial stress (Fig. 3). In this temperature range, the tensile stress accumulated (12–18 MPa) induced failure, resulting in the creation of a through-going fracture and accompanying stress drop. With further cooling, the fracture widened.

Fig. 2 Columnar-jointed basalt at Seljavellir. a Columnar joint outcrop locality complemented by a close-up photo and sketch of striae along a colonnade. The exposure is characterised by quadratic to heptagonal cross-sectional patterns. The fracture surfaces reveal striae, exhibiting both a rough and a smooth portion. b Geometrical relationship between the height of Striae (h) and the column width (d cold ) as shown schematically in Fig. 1. The data are plotted against other columnar-jointed lavas from Columbia River16,17, Staffa5, First Wachung, Prehistoric Mahipuhki and Boiling Pots24 Full size image

Fig. 3 Mechanics of columnar jointing. Tensile stress builds up as a rock, locked into a fixed length, cools at a set rate of 0.05 (orange), 0.1 (red), 1 (black) and 10 °C min−1 (blue), starting with a no applied normal stress and b 5 MPa normal stress. The dashed blue line denotes the solidus temperature (980 °C) Full size image

Dilatometric measurements and mechanical testing

To ensure that this temperature window of columnar jointing is realistic, we support our analysis with a dilatometric and mechanical study to assess whether the dynamics of columnar jointing can be explained by comparison of two distinct, static tests. Dilatometric measurements revealed that the expansion coefficient, α, of the basalt tested is isotropic and linear in the temperature range of interest (400–980 °C), equating to ~10−5 °C−1, with rapid volumetric expansion at temperatures above 980 °C, which indicates melting and matches the solidus temperature estimated by MELTS30 (Fig. 4). Additionally, ambient and high-temperature (820–900 °C) uniaxial compressive strength tests conducted at a strain rate of 10−5 s−1 were used to define the temperature-dependence of the Young’s modulus E (Fig. 5a). E was found to evolve according to an empirical relationship \({\mathbf{E}} = p{\mathbf{T}} + E_0\) (where E is in Pa and T is in °C). Here p = 1.371 × 107 °C−1 and E 0 = 2.6157 × 1010 Pa in the high-temperature window of columnar jointing (Fig. 5b). Finally, ambient and high-temperature (600–940 °C) Brazilian tensile tests constrained the tensile strength of our samples to 12–21 MPa (Supplementary Fig. 4), in good agreement with the failure stress of 12–18 MPa recorded in the combined thermo-mechanical jointing experiments. Together, these thermo-mechanical constraints allow us to model the tensile stress build-up over a range of undercooling ΔT, via \({\mathbf{\sigma }}_{\mathbf{T}} = {\mathbf{E}}\alpha \Delta T\)31. Our calculations suggest that the temperature of macroscopic failure, T f , is 87–144 °C below the solidus; that is, between 893 and 836 °C (Fig. 6)—a temperature window in excellent agreement with the results of columnar jointing experiments developed here (890–840 °C; Fig. 3).

Fig. 4 Thermal expansion coefficient. Dilatometric measurements showing the linear expansion coefficient of Seljavellir basalt, cored axially and parallel to the column, during cooling from 1020 °C. The dashed blue line denotes the solidus temperature (980 °C) Full size image

Fig. 5 Young’s modulus variation with temperature. a Mechanical data obtained through uniaxial compressive testing at ambient temperature (black), 820 (blue), 870 (magenta) and 900 °C (red). Young’s moduli at the different temperatures (E 25 , E 820 , E 870 and E 900 ) are calculated using the linear (elastic) portion of each curve. b Young’s modulus values at each temperature, showing a linear evolution at high temperatures, as depicted by the dotted black line. The dashed blue line denotes the solidus temperature (980 °C) Full size image

Fig. 6 Predicted columnar jointing temperature window. The horizontal dotted black lines show the minimum and maximum strength limits obtained through Brazilian tensile tests. The solid blue line shows the calculated tensile stress build-up upon cooling. The intercepts (vertical dotted lines) between the calculated stress curve and the measured strength lines denote the cooling window necessary to achieve failure. The coloured dots represent the fracturing temperatures achieved during columnar jointing experiments. The dashed blue line denotes the solidus temperature (980 °C) Full size image

Modelling of fracture widening and fluid flow

Studies assessing the thermal history of magma reservoirs, sills and dykes often point towards arguably long cooling timescales if conduction alone is considered, and thus commonly infer the need for external fluid infiltration to increase the cooling efficiency of the magma body2,7,8. The thermo-mechanical constraints introduced here suggest that fluid infiltration may contribute to the thermal budget following columnar joint formation. Here the data reveal that after their formation at 890–840 °C, fractures would open by continued contraction proportional to an expansion coefficient of ~10−5 °C−1 down to 400 °C (and at a slower rate below this temperature; Fig. 3), thereby constructing the permeable network that allows fluid infiltration in an otherwise largely impermeable rock (measured at 5 × 10−20 m2 without fractures). Using the linear expansion coefficient α established for different temperatures (T) by dilatometric measurements, we model the evolution of fracture width, w due to contraction of the rock (between the centre of two columns with length d hot ) when cooling below the fracturing temperature T f (893–836 °C):

$$w\, = \,\alpha d_{\mathrm{hot}}\left[ {T_{\mathrm{f}}\, - \,T} \right]$$ (1)

Here, we gauge the spectrum of d hot using the range of column sizes observed at Seljavellir (d cold = 0.3–1.3 m) as proxy. Our calculation of the final fracture width shows good agreement with values of 1.9 and 8.3 mm measured in the field (Fig. 7a). Knowing the fracture width at different temperatures allows us to predict the permeability of a fractured rock mass, κ fr (in m2), contracting during cooling in the absence of stress, using a scaling for the permeability of fractured systems32. We modify this scaling to account for arrays of hexagonal columnar joints separated by fractures (see Methods) to find:

$$\kappa _{\mathrm{fr}}\, = \,\frac{{s\sqrt 3 \kappa _{\mathrm{h}}}}{{2w + s\sqrt 3 }}\, + \,\frac{{w^3}}{{12w + 6s\sqrt 3 }}$$ (2)

where κ h is the permeability of the intact rock before jointing and s is the edge-length of the colonnade hexagons (in m). Thus, we calculate the permeability of columnar-jointed magmatic aureoles during cooling and show that the evolution is primarily, and highly dependent on the column size (Fig. 7b). Hence the formation of columnar joints may strongly influence hydrothermal circulation and therefore, the cooling history of magmatic bodies near the Earth’s surface and in fluid-rich environments characteristic of shallow volcanic, geothermal and hydrothermal systems.