The Permeability of Columnar Jointed Lava

Columnar jointed lava is an important facies in many geothermal reservoir systems. The permeability of jointed lavas is dominated by the contribution from fracture networks. We use a scaling for the permeability of a set of fractures in a solid or porous mass and extend this to arrays of hexagonal intercolumn fractures. To validate our analytical results, we create numerical domains with relevant geometries and extract system‐scale permeability using the LBflow lattice‐Boltzmann fluid flow simulation tool. Finally, we model the cooling contraction of columns to extend our results so that they predict the permeability with time after lava emplacement. Importantly, we use these results to estimate the range of permeabilities typical of columnar joints that form during cooling from high temperature and are preserved in the crust at moderate to low temperatures.


Introduction
Columnar joints are spectacular features of large mafic lavas, and the dynamics of their formation or their structure and morphology have been the subject of a large body of research to date (Budkewitsch & Robin, 1994;DeGraff & Aydin, 1987;Goehring, 2013;Goehring et al., 2006;Goehring & Morris, 2005, 2008Hetényi et al., 2012;Lamur et al., 2018;Long & Wood, 1986;Müller, 1998bMüller, , 1998aPeck & Minakami, 1968;Phillips et al., 2013). Columnar jointed lavas can be important geothermal reservoir rocks (Fridleifsson & Elders, 2005;Hardardóttir et al., 2009). At present, it is not clear if these jointed basalt lavas are net barriers to flow or effective reservoirs hosting prodigious fluid flow. The Darcian permeability of the bulk of the solid rock samples collected from such lavas is measurable easily in laboratory permeameters and is generally low compared with sedimentary rocks (Lamur et al., 2018;Saar & Manga, 1999). However, the fractures formed by jointing are often wide and open-a result of both cooling contraction after column formation, and erosion during weathering-and not amenable to direct fluid-flow measurements from which permeability can be obtained. In such scenarios, scaling for the permeability of large fractures is needed (Farquharson & Wadsworth, 2018;Heap & Kennedy, 2016;Zimmerman & Bodvarsson, 1996). Here, we develop a framework for the permeability of hexagonally jointed lavas as they cool.
Joints in basalt lavas can occur in a variety of geometries (Goehring, 2013;Long & Wood, 1986;Müller, 1998bMüller, , 1998aPhillips et al., 2013), with the number of sides per column varying between three and eight (Müller, 1998b;Phillips et al., 2013;Walker, 1993). In Figure 1, we compile data for the distribution of sides in natural basalt columns (Budkewitsch & Robin, 1994;Hetényi et al., 2012;Phillips et al., 2013), numerical simulations (Budkewitsch & Robin, 1994;Crain, 1978;Tanemura & Hasegawa, 1980), and analog experiments in drying corn starch systems (Müller, 1998b(Müller, , 1998a, confirming that within this variability, six-sided columns are the most frequent. As well as being six-sided, for our purposes, it is important to know the regularity of the polygonal shapes. To test this, we compile data for the side-length of the six-sided shape b, the two-dimensional diameter of the polygonal top D, and the hexagonal area (Hetényi et al., 2012;Phillips et al., 2013). If the columns are regular hexagons, these parameters should scale as D ¼ ffiffi ffi 3 p b, which appears to be robust across a wide range of scales (Figure 1).
In the colonnade of columnar jointed lavas, the columnar nature of the jointing also appears to be restricted to a two-dimensional plane, in which the long-axis of the fracture arrays are all parallel, broadly perpendicular to a lava top and bottom surfaces (Long & Wood, 1986;Phillips et al., 2013;Walker, 1993). For the reasons given here (Figure 1), in what follows, we concentrate our analysis on regular hexagonal geometries in two dimensions. Based on the success we find with this geometry, we provide the equivalent scaling for unit cells of other tessellating polygonal arrays and explore the evolution of permeability during cooling.

Scaling Approach for the Permeability of Fractured Lava
In this section, we provide a derived analytical scaling law for the permeability as a function of the porosity of basaltic lavas fractured in perfect hexagonal arrays. We start with the permeability of a single fracture in an infinite medium and extend that result to the permeability of many fractures in a system of any size.

The Permeability in the Direction Parallel a Single Smooth Fracture in an Infinite Medium
First, we must find the permeability of a single smooth fracture of infinite length but finite width, w. To do this, we use the derived form of Poiseuille's law for the steady-state velocity distribution of a fluid in a slot flowing down a pressure gradient ∇P in the low-Reynolds number regime, when no-slip conditions exist at the walls of slot, which has the form where u(x) is the velocity in the axial direction, which varies with x-position (orthogonal direction), and μ is the fluid viscosity. Integrating u(x) to find the gap-averaged velocity 〈u〉, yields (Langlois & Deville, 2014;Witherspoon et al., 1980) where we use 〈·〉 to denote any average property across the gap. Using Darcy's law to relate the general fluid velocity tensor 〈u〉 to a permeability tensor k, of the form 〈u〉 = − k ∇ P/μ, we can find the generic solution for the permeability in the axial direction and in smooth slots, termed k ′, as which provides the convenient dimensional scaling for fracture permeability in the axial direction, as k≈k 0 =w 2 , where a bar above a parameter denotes it is dimensionless. This finding is consistent with previous work on the permeability in fractured systems (Heap & Kennedy, 2016;Zimmerman & Bodvarsson, 1996) but given here for completeness.  Hetényi et al. (2012), Müller (1998b), and Budkewitsch and Robin (1994) and references therein. (c) The relationship between column diameter D and the side length of a column b, for hexagonal columns (six sides) scaling as D ¼ ffiffi ffi 3 p b. Inset-the relationship between column area and side length for hexagonal columns (six sides) scaling as A ¼ 3 ffiffi ffi 3 p b 2 =2 and pentagonal columns (five sides) scaling as The data in panel (c) are from published sources (Hetényi et al., 2012;Müller, 1998aMüller, , 1998bPeck & Minakami, 1968;

The Permeability in the Direction Parallel to an Array of Smooth Fractures
The result given in equation (3) is not directly useful for determining the permeability of any real system as it requires scaling against a system-size. Using 〈k〉 to denote the permeability of the system, classic results for this scaling are for the case where fractures are parallel to one another and of equal width can be scaled by 〈k〉 = ϕ f k where ϕ f = A f /A T is fracture surface fraction (i.e., the ratio of the fracture surface area to the total surface area of the system; Lucia, 1983). The surface areas are measured in 2D, rather than the direct surfaces of all the fractures. This assumes the nonfractured host rock is impermeable to fluid. If we relax that assumption, we can add the scaled contribution from the host by (Ebigbo et al., 2016;Heap & Kennedy, 2016) where k f is the total fracture permeability, ϕ h is the host surface fraction, k h is the host permeability, which we take to be isotropic, and A h is the cross-sectional area of the host rock, such that A h = A T − A f . When the fractures are of nonequal width, the contributions from each width class can be summed and weighted by the area of that class, following solutions for directionally dependent permeability in anisotropic systems (Farquharson & Wadsworth, 2018). We note that this is different from the approach of Lamur et al. (2017) in which k h is added to k′ without scaling the contributions. In what follows, we develop equation (4) for the specific case of hexagonal columns formed in cooled lava.

The Permeability of Hexagonally Jointed Lava
Now we can turn to the system at hand, which is, in essence, hexagonal, for the reasons given ( Figure 1) and following evidence from natural basaltic columnar joint geometries (Budkewitsch & Robin, 1994;Hetényi et al., 2012;Müller, 1998b;Phillips et al., 2013;Walker, 1993). Now A T is the sum of the areas of a set of tessellating hexagons and the areas of all the fractures. If we assume that the sides of the hexagons are long compared with the widths of the fractures, which is reasonable for columnar joints, then we can, to a first approximation, ignore the complications of the fracture vertices and approximate the fracture area A f as being a sum of rectangular fractures (see our tessellating grid in Figure 2a).

The area of each of the hexagons is
where b is the side length of one hexagon. Then if we take a unit cell of a tessellating hexagon array (Figure 2a), we have that A h is that of a single hexagonal area, so The unit cell contains three fracture rectangles each of which with an area of The region of interest used with LBflow is smaller than the total domain and, for each repeat measurement, is moved around the larger simulation grid. Here we have oriented the domain such that the "flat topped" orientation is shown (imagine rotating the image 90°and viewing the "pointy topped" hexagonal array). Inset-a "pointy topped" zoom-in of a portion of the total array, showing the "unit cell" used for the scaling of hexagonal properties. (b) Calibrating the critical size of the region of interest, in order to achieve reproducible and accurate results. Here plotted is the output permeability as a function of the domain volume L 3 showing that domain sizes ≳10 −12 m two possibilities: in the x-direction, the hexagons are arranged in a "pointy topped" array and in the y-direction in a "flat topped" array. We can use the rotation matrix method (described in the context of volcanic rocks in Farquharson & Wadsworth, 2018) to scale the two-dimensional tensorial description of the permeability of a hexagonal array in each mutually perpendicular horizontal direction, which we term b k, as follows where k is the permeability tensor prior to rotation and R is the rotation matrix, such that For an unrotated fracture whose long axis is parallel to the x-axis, one can set k 1 = k' and k 2 = 0, which implies that equation (5) and so b k x ¼ k ′ cos 2 θ and b k y ¼ k ′ sin 2 θ. We can now add the scaled in-series contributions of each component of a set of fractures in hexagonal arrays via k ′ sin 2 θ n where n represents the nth fracture, N the total number of fractures and θ n the angle of the nth fracture with respect to the x-axis. In the unit cell, one fracture rectangle is parallel to the x-axis (θ = 0) while the other two are rotated by an angle of θ = π/3 and θ = − π/3.
. Therefore, in both the x-and ydirections, we have the same scaling, even though the spatial distribution of local flow rates will be different.
For flow in the z-direction, the fracture contribution is simply A f k ' /A T . Each fracture has a permeability k ′ = w 2 /12. This allows us to rewrite equation (4) as an analytical expression for the permeability of an array of columnar joints with ignoring the contributions of the corner segments of the fracture network, which are triangular in cross section ( Figure 2a).
As a complicating step, we can take account of the triangular vertices, two of which appear in each unit cell.
For flow in the z-direction equation (4) becomes is the area of the vertices in the unit cell (which has to be incorporated in A T = A f +A v +A h as well when solving equation (9)), and k v = w 2 /80 is the permeability of a triangular duct (adapted from the conductance; Llewellin, 2010b;Sisavath et al., 2001). This expression for k v is true for the case here where the side length of the triangular vertex is equal to the fracture widths of the rectangular fractures. Strictly speaking, this is not valid, because the walls of the triangular duct are not solid and are, instead, the terminating ends of the rectangular portions of the fractures. However, we include this as an approximation and is very much a lower limit for the effect of the vertices. For brevity, we omit to write out the solution to equation (9) in full equivalent form as given in equation (8) without vertices.
Equation (9) is given as a function of the ratio of two-dimensional cross-sectional areas. Another form is to give it as a function of volumetric fractions. This is done by noting that the third dimension lengthscale (the depth of the fractures) is the same for the fractures, vertices, and host rock, and therefore, the volumetric version of equation (9) is simply where ϕ v = A v /A T and ϕ h is the total volume fraction of the host rock (and not the porosity of the host rock). Under the assumption that the host rock is negligibly porous, the total pore volume fraction of the system would be ϕ T = ϕ f +ϕ v for equation (10) and ϕ T = ϕ f for equation (8). We note that this assumption of a nonporous host rock is not required and that equations (4), (8), (9), and (10) are valid for any case.

Numerical Methods
The analytical solutions found for the permeability in equations (8) and (10) are appealing in their simplicity. However, some assumptions were made in the derivation that require validation before applying them to large systems. To validate the solution, one could perform carefully constructed experiments or numerical tests. Given the challenges associated with making precise geometries at small sample scales and in maintaining through-going fractures at defined widths, we opt for the latter approach. Using Python™, we create an algorithm that populates a numerical cubic box (edge length L) with a hexagonal array in the geometry given in Figure 2a. We set up our simulation domain such that the center of each hexagon is defined, meaning that as we control the edge length b, w evolves proportionally. Our control parameter is therefore w.
For each domain that we create in three dimensions, we define the fractures as being fluid-filled, and the columns as being solid, and then we use the LBflow (Llewellin, 2010a(Llewellin, , 2010b code to simulate a fluid flow through the fractures. The LBflow code is versatile and has been validated against analytical results for the flow between lattices of nonoverlapping spheres in body-centered-cubic (Llewellin, 2010b), and both simplecubic and face-centered-cubic configurations . It has also been used as an exploratory tool for determining the permeability between fibrous media (Nabovati et al., 2009) and between particles in dynamic sintering experiments . The code discretizes the fluid phase into a cubic lattice of nodes (D 3 Q 15 lattice), and packets of mass are propagated with time t (see Llewellin, 2010a). We use a simulation fluid with the properties μ = 1.8205 × 10 −5 Pa s, and density ρ = 1.2047 kg m 3 . The simulation is run for a time sufficient to ensure that steady state has been reached. Steady state is assessed as being the time when the magnitude of the average velocity of the fluid nodes u converges to within 10 −5 of the mean value for 50 time steps two times consecutively. The magnitude of u is then converted to a volume average by 〈u〉 = uϕ, where ϕ is the volume fraction of the fluid in the total system. The conditions are checked such that they must satisfy low-Reynolds number creeping flow and low Mach number flow. Defining these as Re = ρL〈u〉/μ and Ma = 〈u〉/c, respectively, where L is the length of the system in the direction of flow, and c is the sound speed in the fluid, we find that we cover the range 10 −10 < Re < 10 −7 and 10 −14 < Ma < 10 −10 . We use an imposed pressure gradient of ∇P = 0.01 Pa.m -1 .
We ran calibration numerical experiments to check for domain size effects. We found that for domain volumes of ≤10 6 px 3 , repeat runs on randomly realized domains were not reproducible. For domain volumes >10 −12 m 3 , repeat runs were reproducible and converged on the theoretical result for flow in the z-direction (given above; Figure 2b). We therefore conducted all final simulations on domains with volumes~10 −11 m 3 to ensure scale-independent results.
Our simulation results output a value of 〈u〉 for each domain (Figure 2c), which, knowing that Re ≪ 1, we can convert to a permeability using Darcy's law, stated earlier. For simplicity, we set the solid columns to be impermeable. In detail, we subsample a larger domain of hexagonal arrays in order to run the fluid-flow simulation. To assure that we are sampling a representative size of the full domain, we run several subsamples per simulation condition, averaging the result and using the standard error on the average as a means to compute the uncertainty via the standard error calculation. The small uncertainty over a large number of repeat measurements on various subsamples gives us confidence that our solution captures the large-scale behavior and is therefore scale-independent. In Figure 3, we show a representative visualization (using Paraview™) of the fluid flow speed distribution in each principle direction.

Results and Validation
The output of our simulations is a measured 〈k〉 in all three principal directions, which can be compared with the results of equations (8) and (10) directly. When comparing the dimensional permeability (in m 2 ), we find excellent agreement between our simulation results and the analytical result given in equations (8) and (10), across a wide range of w (Figure 4a). We note that for the simulation results, k h = 0, meaning that the validation is strictly for the reduced form of equations (8) and (10), 〈k x 〉 = 〈k y 〉 = ϕ f k ' /2 = w 3 /(24w +12 √3b) and 〈k z 〉 = ϕ f k′ = w 3 /(12w+6 √3b), and 〈k z 〉 = ϕ f k′ +ϕ v k v , respectively.
Similarly, if we scale 〈k〉 by the equivalent fracture permeability k e = k′ and k e = (ϕ f k′+ϕ v k v )/(ϕ f +ϕ v ) in equations (8) and (10) respectively, which yields k k e ¼ ϕ T , and if we replace w with the gas volume fraction in all the fractures in the system ϕ T , we find a universal dimensionless solution, which we propose is valid for predicting 〈k〉 across a wide range of w and b value pairs (Figure 4b). We provide all simulation results in Table 1.

Discussion
Here we discuss a few extensions to the analysis given above, which may be particularly useful for interpreting fluid flow in fractured lava systems that have a column-like geometry, meaning that they can be assumed to have all fractures parallel to one another. For these extensions, we do not explicitly validate the results, relying instead on the efficacy of the main analysis given above, as justification for these ideas.

Flow Through Nonhexagonally Jointed Lava
Basaltic columns are most commonly jointed in hexagonal arrays of fractures, particularly in the colonnade section(s) (Müller, 1998b;Phillips et al., 2013;Walker, 1993). However, there exist a variety of nonhexagonal sections. In this paper, we have validated a scaling law for hexagonal columns, using numerical methods. This gives us some confidence that similar assumptions can be made to arrive at scaling laws for other jointing patterns polygonal in two-dimensional cross section.
First, we can consider the three possible monohedral tessellating regular polygon tiles in two-dimensional cross section, which include not only the hexagonal tile but also the equilateral triangle and square. Following the approach given in previous sections, for the latter two polygons, assuming that in both cases, the regular polygon edge length is always termed b, we have variations of equation (8) as for a triangular tile, and for a square tile, where the approximation of the fractures as sets of angled cuboids all with width w leaves out hexagonal and square shaped vertices for the triangular and square tiles, respectively, which are not accounted for by equations (11) and (12).

The Permeability During Cooling Contraction of Fractured Lava
Returning to the hexagonal geometry, it is clear from previous work on the formation and growth of basaltic columnar joints that these fractures form while hot and grow during cooling (Budkewitsch & Robin, 1994;Lamur et al., 2018). Therefore, both b and w are functions of time during cooling. In what follows, we do not address the fracture propagation and growth dynamics and instead take a simplified approach to predict w(t) for a range of initial column edge lengths at the time of formation b i .
We assume that the columns are purely elastic solids, which is implied by recent work (Lamur et al., 2018). In this case, the area of the columns will reduce as temperature reduces by an amount ΔT, opening up the fractures between them as (Turcotte & Schubert, 2014) where A h,i is the initial area of the column at the point of formation, and α A is the area expansivity (in K −1 ). Equation (13) is only strictly true if the expansivity is broadly constant with temperature and if the expansion can be assumed to be linear. Area expansivities are rarely measured directly but can be related to the more commonly measured linear expansivity α L for isotropic solids by α A ≈ 2α L . The range of values for α L Figure 4. (a) The results for the permeability 〈k x 〉, 〈k y 〉,and 〈k z 〉 of a hexagonally fractured system as a function of ϕ T . (b) The same results as in a, but for the normalized k . In both panels, the curves represent the solutions of our analytical scaling (solid is for hexagonally arranged rectangular fractures of width w using equation (8), and dashed is the same but with the addition of the triangular vertices using equation (10)). When the results are rendered dimensionless in b, there is a single solution for k z (top curve) and a single solution for k x ¼ k y (bottom curve). Inset-the same data, but where ϕ T is converted to w/b. Each data point shape is for a different hexagon side length b, marked here as b i because in detail, the simulations were performed by decreasing the hexagon side length in order to increase ϕ T , rather than by directly controlling w; for which the results would be the same. measured for basalt is approximately 1 × 10 −5 K −1 between 200 < T < 980°C (Lamur et al., 2018). For simplicity, we assume this is extensible down to ambient conditions T~20°C. Furthermore, Lamur et al. (2018) show that the minimum temperature below, which contraction occurs after fracture formation, can be taken as 840°C. This gives a significant fractional area change of ΔA h /A h,i = 0.0164 for ΔT = 820°C.
By way of a worked example, we can convert ΔA h /A h,i to 〈k z 〉 using the tools developed herein. First, we note that there are abundant data for the geometry of the columns (e.g., b and A h ) but scarce data for the fracture widths w between columns. Therefore, it is useful to find a solution for w as a function of the final geometry we observe. We find that using equation (13) and the geometry of hexagonal columns, we can relate the hexagonal edge length measured in the field after the rock has cooled, termed b f , to the initial edge length at the temperature at which the columns formed, termed b i , by b i = b f (1+α A ΔT) 1/2 . In turn, this allows us to compute the value of w that results from the cooling over the interval ΔT, by w ¼ ffiffi ffi Using the observation that 0.1 < b f < 1 m (Figure 1) and that a typical ΔT is 820°C to find that an approximate solution is 0.00144 < w < 0.0144 m, and therefore, using equation (8), we have that 2.83 × 10 −9 < 〈k z 〉 < 2.83 × 10 −7 m 2 . We propose that this is a reasonable scaled estimate for columnar jointed lava that could be used in the absence of more information.

Inertial Contributions to Fluid Flow
In the above, we have limited ourselves to Darcian, low-Reynolds number flow, and the Darcian permeability that can be determined as a result. However, in the crust, low viscosity fluids and gases have the potential for inertial contributions to flow. In this case, the bulk flow velocity is where k I is the inertial permeability tensor. The LBflow tool is not specifically validated for turbulent flow at high-Reynolds number, and therefore our simulation results were specifically for the low-Reynolds number end member case. However, Zhou et al. (2019) have demonstrated that across a very wide range of materials of geological relevance, including fractured structures such as are relevant here, k I ∝ k a , where a is a constant. More specifically, they propose with ω ≈ 10 10 m 1−α and α ≈ 3. In this study, we validate a law for the full tensor of k in columnar jointed lavas, and we propose that equation (15) provides the tool required to unify this result across all Reynolds numbers in the low-Mach number regime.
To give an example of the utility of our approach to constraining the permeability (equation (8)) and the inertial permeability (equation (15)), we present a solution to the full Forchheimer equation (equation (14)) for water flowing through hexagonal fractures perpendicular to the face orientation ( Figure 5). This solution shows the evolution of the flow velocity with increasing driving pressure gradient for a range of relevant fracture widths. The transition from Darcian (viscous) flow to inertial flow occurs at a Forchheimer number Fo (which is a modified Reynolds number). To give an example in one direction, the Forchheimer number for flow in the z-direction would be As a first-order approximation, the critical Forchheimer number for the transition to inertial flow would occur around Fo ≈ 1, such that at Fo ≪ 1, flow would be viscous and Darcian, and at Fo ≫ 1 flow would be inertial, turbulent, and non-Darcian. In Figure 5, we give this transition point, showing that it is different for each fracture width.

Limitations and Conclusions
Important limitations of our approach exist, which are worth stating. First, our idealized numerical set up is necessarily constrained to smooth-edged fractures and to cases where the columns are straight in the z-direction throughout the columnar jointed lava. Neither of these are necessarily true in nature. Observations of curved column faces, especially in entablature zones (Hetényi et al., 2012), are necessarily neglected here. We anticipate that fracture surface roughness is important when the fractures are narrow (Heap & Kennedy, 2016) and can be neglected in very wide fractures so long as the roughness lengthscale is much less than the fracture width (Brush & Thomson, 2003;Mallikamas & Rajaram, 2010;Wang et al., 2015). This is consistent with the result of the large datasets compiled in Zhou et al. (2019), wherein it is shown that surface roughness exerts a seemingly second-order effect.
We present a physical scale for predicting the permeability of jointed lavas. We explore this explicitly for hexagonally jointed systems using numerical methods and a lattice-Boltzmann fluid flow simulation. Finally, we use the efficacy of our scaling to explore alternative and more general geometries for lavas that might be fractured in less regular arrays. Our result suggests that images of the top surface of cold lavas could be easily analyzed to extract the area fraction of the fractures and their distribution of widths and lengths in order to predict the permeability of the entire rock mass for the case when the pressure potential driving fluid flow is in the direction of the long-axis of the fracture network. This approach is useful for geothermal systems that comprise lava as a significant lithofacies. Figure 5. The flow velocity in the z-direction (perpendicular to the polygonal face orientation) 〈u z 〉, in a hexagonally jointed lava as a function of the fluid pressure gradient ∇P using equation (14) (with equation (15) to predict k I ). Here the pore fluid is assumed to be water with properties μ ≈ 10 −3 Pa. s and ρ = 1000 kg. m −1 . The host rock permeability is set to zero k h = 0 (valid for nominally impermeable basalt; c.f. measurements k h ≈ 10 −20 m 2 ; Lamur et al., 2018). The side length of the hexagons modelled here is b = 1 m.We give the result for fracture widths w = 10 −3 m and w = 10 −1 m as dashed lines, as well as the results for every width between these end members as a continuous color mapping (scale indicated on the plot). The dotted line represents the approximate transition from Darcian (or viscous) flow to inertial flow at Fo = 1 (see text for details).