Amalie L Christensen

Christophe Raufaste

Marek Misztal

Franck Celestini

Maria Guidi

Clive Ellegaard

Joachim Mathiesen

Scale selection in columnar jointing: insights from experiments on cooling stearic acid and numerical simulations

Scale selection in columnar jointing is a nontrivial function of material properties and cooling conditions • Columnar jointing is reproduced experimentally in a model system using cooling stearic acid • Proposed model can be used in the prediction of the cooling rate of columnar jointed systems

Introduction

Spectacular geometric patterns, such as Giant's Causeway at the Irish shoreline, Devil's Tower in Wyoming or Svartifoss in Iceland (Figure 1a), are examples of columnar jointing in igneous rocks. Several hypotheses for the physical process resulting in columnar jointing have historically been proposed, such as the columns being crystals of quartz or the columns being regular packings of spherulites [START_REF] Tomkeieff | The basalt lavas of the Giant's Causeway district of Northern Ireland[END_REF]. It is now widely accepted that columnar jointing in igneous rock is a result of thermal contraction, an idea developed in the late 19th century [START_REF] Mallet | On the origin and mechanism of production of the prismatic (or columnar) structure of basalt[END_REF]. As the igneous body cools it contracts. If the contraction is uniform, no stress will be generated in the rock, whereas a non-uniform contraction, e.g. due to temperature gradients, will generate stress. Above the glass transition temperature of the rock, T G , the stress is viscously dissipated, but below T G stress starts to accumulate. When the stress exceeds the fracture threshold, the accumulated elastic energy is transformed to surface energy by fractures. As the rock further cools with time, the T G isotherm propagates toward the interior of the rock and with it the fracture front. In that way a network of fractures, or so-called joints, is created and it divides the rock into columns.

Igneous columnar jointing occurs in lava flows, lava lakes, sills, dikes and ash-flow tuffs [START_REF] Hetényi | Scales of columnar jointing in igneous rocks: field measurements and controlling factors[END_REF][START_REF] Spry | The origin of columnar jointing, particularly in basalt flows[END_REF]. The columns are commonly polygonal shaped, with the diameter being approximately constant over a significant fraction of their height. The observed diameters range from a few cm up to 3 m while the column height can be up to 30 m. The regular columns tend to have 5, 6 or 7 sides, although 3-, 4-, and 8-sided columns exist as well [START_REF] Hetényi | Scales of columnar jointing in igneous rocks: field measurements and controlling factors[END_REF][START_REF] Phillips | The formation of columnar joints produced by cooling in basalt at Staffa, Scotland[END_REF][START_REF] Spry | The origin of columnar jointing, particularly in basalt flows[END_REF]. As illustrated D R A F T February 6, 2022, 10:10am D R A F T in Figure 1b, the regular columns can be intersected by a highly disorganized region with smaller curvy columns. The two types of regions are commonly termed colonnade and entablature, respectively [START_REF] Tomkeieff | The basalt lavas of the Giant's Causeway district of Northern Ireland[END_REF]. Figure 1c illustrates striations on the surface of columnar joints. The striations are linear bands oriented perpendicular to the crack propagation direction and are produced by stepwise crack advance, where the crack advance length corresponds to the striae height, s [START_REF] Ryan | The glass transition in basalt[END_REF].

Under certain conditions, columnar jointing can be reproduced at the lab scale in desiccating starch slurries [START_REF] Goehring | Desiccation Cracks and their Patterns: Formation and Modelling in Science and Nature[END_REF][START_REF] Goehring | Experimental investigation of the scaling of columnar joints[END_REF][START_REF] Müller | Experimental simulation of basalt columns[END_REF][START_REF] Müller | Experimental simulation of joint morphology[END_REF][START_REF] Toramaru | Columnar joint morphology and cooling rate: A starch-water mixture experiment[END_REF]. In general, desiccation cracks differ from columnar jointing, since the contraction of the desiccating substance, e.g. mud, is limited to the surface in contact with air. The cracks therefore have a short penetration depth [START_REF] Groisman | An Experimental Study of Cracking Induced by Desiccation[END_REF]. In contrast to columnar jointing, mud cracks typically have a hierarchical structure, where domains bounded by existing cracks are divided by secondary cracks meeting existing cracks at 90 degree angles [START_REF] Bohn | Four Sided Domains in Hierarchical Space Dividing Patterns[END_REF]. Hierarchical fracture patterns are also produced by serpentinization processes in e.g. dyke-host complexes [START_REF] Iyer | Reactionassisted hierarchical fracturing during serpentinization[END_REF] and by weathering processes [START_REF] Røyne | Controls on rock weathering rates by reaction-induced hierarchical fracturing[END_REF]. Common for most hierarchical fracture patterns is the lack of a maturation process similar to that observed in cooling basalt, where the fracture network has time to mature as it propagates with the cooling front [START_REF] Hofmann | Why Hexagonal Basalt Columns?[END_REF].

The internal stress generated by the contraction of a material, β = α T ∆T , will in a linear elastic approximation be proportional to Eβ/(1 -2ν), where E is the Young modulus, α T the coefficient of linear thermal expansion, ∆T a temperature change and ν the Poisson ratio. The stress generation will be lower if the material contraction changes from 0 to D R A F T February 6, 2022, 10:10am D R A F T β over a length scale w and in the limit of large w, the stress becomes vanishingly small.

Inspired by desiccating thin films and mud, we write the magnitude of the internal stress σ in terms of a scaling function f :

σ = Eβf (w, ν, . . .) (1)
The internal stress is dependent on material properties, such as Young's modulus E,

Poisson's ratio ν and coefficient of linear thermal expansion α T , as well as on the externally imposed cooling conditions determined by the temperature change ∆T and the length scale w. The internal stress might depend on further material properties such as the thermal diffusivity D, the fracture toughness K I,c or cooling conditions.

Asymptotic shapes of such a scaling function has been suggested in [START_REF] Hayakawa | Pattern selection of multicrack propagation in quenched crystals[END_REF].

For desiccation in thin film systems, the scaling function can be estimated analytically [START_REF] Cohen | Drying patterns: Sensitivity to residual stresses[END_REF]. Here, we will determine the shape and dependencies of the scaling function, f , for columnar jointing from simple analytical models and from numerical simulations and compare with experiments on cooling stearic acid.

Heat extraction mechanism

The stress accumulation in a cooling igneous body is controlled by the magnitude of the temperature gradients, which again is controlled by the heat transport in the material. Two main modes of heat transport during columnar jointing have been considered

in the literature, bulk heat conduction and crack aided convective cooling, respectively.

In general, conductive cooling leads to a temperature profile flattening as time progresses and therefore results in an increase of column diameters with time. In contrast, convective cooling can maintain a steep temperature profile and therefore also a fixed column

D R A F T February 6, 2022, 10:10am D R A F T
diameter. Below we shall focus on the convective cooling mode and only briefly discuss the conductive mode.

Conductive cooling

As a simple model, one can consider the igneous body as a homogeneous and isotropic half space X ≥ 0 where we imagine the X-axis to point into the igneous body. Due to the assumed homogeneity and infinite extent of the Y -and Z-dimensions, the temperature field is only a function of X. For pure bulk heat diffusion the temperature field is governed by:

∂T (X, t) ∂t = D∇ 2 T (X, t) (2)
At the interface, X = 0, the igneous body is in contact with a cooling surface which we assume efficiently removes heat, such that the temperature of this surface is always T 0 , i.e. T (X = 0, t) = T 0 . The initial temperature of the body at t = 0 is taken to be homogeneously T 1 , i.e. T (X > 0, t = 0) = T 1 . If we do not consider latent heat, the present boundary conditions yield:

T (X, t) = T 0 + (T 1 -T 0 ) erf X 2 √ Dt (3)
where the error function is defined as: erf(η) = 2 √ π η 0 e -τ 2 dτ . If we assume that the igneous body can be described as a linear elastic material (infinite in the Y -and Zdimensions), then from Eq. ( 3), we see that the only natural length scale in the system must be proportional to √ Dt. We would therefore expect the stria height as well as the average column diameter to increase with the square root of time, or equivalently, the distance from the cooling surface. Here, we have neglected the latent heat, but the result from including latent heat is qualitatively identical. The profile again flattens proportional to the square root of time. While an increase in striae height has been observed within a D R A F T February 6, 2022, 10:10am D R A F T few meters from the cooling surfaces [START_REF] D R A F T Ryan | Cyclic fracture mechanisms in cooling basalt[END_REF][START_REF] Degraff | Effect of thermal regime on growth increment and spacing of contraction joints in basaltic lava[END_REF],

jointed igneous rock commonly has a constant striae height over much of the column height [START_REF] Goehring | Scaling of columnar joints in basalt[END_REF][START_REF] Phillips | The formation of columnar joints produced by cooling in basalt at Staffa, Scotland[END_REF]. The constant stria height could be explained if the system had one fixed length scale. A fixed length scale is a possibility in convective cooled systems.

Convective cooling

Borehole temperature measurements of the Kilauea Iki lava lake [START_REF] Hardee | Solidification in Kilauea Iki lava lake[END_REF] showed a constant temperature of 100 • C in the upper part of the solidified basalt. This suggests, that the mechanism responsible for stabilizing the heat extraction rate and thereby the speed of the temperature profile is convection of water/steam in the joints. We shall therefore assume that the water and steam zone propagates through the cooling body at a constant velocity v = v X and that the entire convection zone is at temperature

T 0 = 100 • C.
The water and steam convection results in a moving boundary condition of T (X, t) = T 0 at the position X -vt = 0 for the heat diffusion equation in Eq. ( 2). To implement this boundary condition, we change the reference to a frame moving with the convection zone and rescale by a characteristic length, , of the system, x = (X -vt)/ .

In steady state, the heat diffusion equation in the moving frame reduces to:

1 P e ∂ 2 T ∂x 2 + ∂T ∂x = 0, (4)
where we have introduced the Péclet number P e = v /D which is the ratio of the heat advection and heat diffusion rates. The Péclet number can also be thought of as the ratio, P e = /w, between the column diameter and the width of the thermal front w = D/v.

D R A F T February 6, 2022, 10:10am D R A F T
The general solution to Eq. ( 4) is:

T (x) = ∆T 1 -e -P e x θ(x) + T 0 , (5)
where θ(x) is the Heaviside step function, ∆T = T 1 -T 0 and T 0 , T 1 are the temperatures at x = 0 and x = ∞ respectively. T 0 is the temperature of the convection zone and T 1 is the emplacement temperature. This temperature profile is similar to the profile employed by [START_REF] Hardee | Solidification in Kilauea Iki lava lake[END_REF] which yielded good agreement with borehole measurements at the Kilauea Iki lava lake. Eq. ( 5) has a fixed shape in the frame of reference moving with the convection zone and is therefore expected to result in columnar joints with constant diameter, , and constant striae height, s.

Scaling of the column diameter

Basalt columns typically have a constant column diameter over most of their height [START_REF] Peck | The Formation of Columnar Joints in the Upper Part of Kilauean Lava Lakes, Hawaii[END_REF][START_REF] Goehring | Scaling of columnar joints in basalt[END_REF][START_REF] Phillips | The formation of columnar joints produced by cooling in basalt at Staffa, Scotland[END_REF]. The speed of the temperature front, which is proportional to the cooling rate, is widely considered the governing quantity of the column diameter [START_REF] Budkewitsch | Modelling the evolution of columnar joints[END_REF][START_REF] Grossenbacher | Conductive cooling of lava: columnar joint diameter and stria width as functions of cooling rate and thermal gradient[END_REF], with faster cooling leading to slenderer columns. As [START_REF] Budkewitsch | Modelling the evolution of columnar joints[END_REF] notes, the rate of heat extraction from a convectional cooled igneous body is equal to the heat transferred from the hot interior to the crack network. inverse of the temperature front speed, 1/v [START_REF] Goehring | Nonequilibrium scale selection mechanism for columnar jointing[END_REF]. From the definition of the Péclet number, it was proposed that igneous columnar jointing proceeds at a constant Péclet number of P e = 0.3 ± 0.2.

Hetényi et al.

[2012] further found, that the geological setting (lava flow, lava lake, lava dome, sill, dyke) exerted control on the columnar diameter, through control of the surfaces where heat can be exchanged with the environment and thereby the cooling rate. When the bodies were geometrically unconstrained (lava flows), they found chemistry to play a role, through its control of material properties and through its control of emplacement geometry, which in turn affects the cooling rate.

Though thermal contraction is widely believed to be the mechanism governing the column formation process, other models where fingering in the solidifying lava governs the column scale, such as large scale constitutional supercooling [START_REF] Guy | Comments on "Basalt columns: Large scale constitutional supercooling? by John Gilman (JVGR, 2009) and presentation of some new data[END_REF][START_REF] Gilman | Basalt columns: Large scale constitutional supercooling?[END_REF] or double-diffusive convection [START_REF] D R A F T Kantha | Basalt fingers' -origin of columnar joints?[END_REF], have been proposed. These models were not found consistent with a small-scale study of a single basalt column [START_REF] Bosshard | Origin of internal flow structures in columnar-jointed basalt from Hrepphólar, Iceland: I. Textural and geochemical characterization[END_REF], and do not fully explain the striking similarity between columnar jointing in basalt and corn starch.

All the above considerations assume a one-to-one relation between material properties, cooling parameters and the resulting column diameter. Alternatively, one might imagine that there is a range of possible column diameters for each set of material properties and cooling parameters. [START_REF] Goehring | Experimental investigation of the scaling of columnar joints[END_REF] found in desiccating corn starch experiments, that for the same initial evaporation rate there existed a column diameter stable over a range of final evaporation rates. [START_REF] Jagla | Sequential fragmentation: The origin of columnar quasihexagonal patterns[END_REF] numerically studied the evolution of a parallel array of cracks in a directionally dried system. They found that the distance D R A F T February 6, 2022, 10:10am D R A F T between the parallel cracks was governed by the typical length scale at the surface where cooling was initiated. In contrast to this finding, [START_REF] Bahr | Diameter of basalt columns derived from fracture mechanics bifurcation analysis[END_REF] deduced from a numerical bifurcation analysis in three dimensions, the existence of a smallest possible column diameter for each value of the crack front propagation velocity.

For columnar jointed igneous rock, there is to our knowledge no clear evidence of a oneto-one relation between material properties (such as E, α T , D, K I,c ), cooling parameters (such as ∆T, w) and the column diameter. Though the studies of column diameter versus geological emplacement and column diameter versus cooling rate show correlations, none of these data sets are strong enough to imply a one-to-one relationship. We will therefore explore the idea of a range of possible columns diameters.

Organization of the article

In this article, we shall investigate the scaling of columnar joints through threedimensional numerical simulations, available field data and a new experimental system based on the cooling of stearic acid. The stearic acid experiments are presented in Section 2. A simple two-dimensional model for columnar jointing is introduced in Section 3, and the idea of a range of possible columns diameters for each system configuration is explored analytically and numerically in Section 4. We use the numerical simulations and analytical arguments to derive a quantitative relation between the Péclet number and material properties and cooling conditions. While the analytical arguments follow from two-dimensional considerations, we show from numerical simulations that the arguments are also valid in three dimensions. Initially, the wax is melted and brought to a temperature of T 1 = 120 • C. The liquid is poured inside a cylindrical aluminum can that is placed inside a coolant liquid, liquid nitrogen at T 0 = -196 • C, and solidification occurs from the can edges towards the bulk of the sample. At the end of the solidification process, the sample is broken into pieces to image the inner structure. A columnar structure appears inside the solid stearic acid at the contact with the aluminum can and extends roughly 2 cm inwards (Figure 2a). On the rest of the sample, no organized structure is observed. At the end of the experiment, the sample exhibits a volume decrease of approximately 10%.

Post-solidified samples are imaged in depth by X-ray microtomography (Figure 2c-2f).

The scanner is a SkyScan1172, with a source voltage of 59 kV and a source current of 167 µA. The stearic acid samples have a surface of order 15 × 15 mm 2 and a depth of about 20 mm. The average area of the columns is analyzed and plotted as a function of the depth from the cooling surface in Figure 2b. We find that the column area is rather constant through the jointed region though an increasing trend might be present.

D R A F T February 6, 2022, 10:10am D R A F T
Assuming perfect hexagonal columns, the average column area is used to estimate the column diameter, = (2 ± 1) mm.

The solidification front velocity is estimated by stopping the solidification of similar samples at different times (liquid stearic acid is removed from the can in less than 2s) and measuring a posteriori the solid thickness. We find that the velocity is rather constant with a typical solidification front velocity of 5 µm/s.

In situ measurements are performed as well by immersing an array of thermocouples (0.4 mm between each of them spanning 5.6 mm in total) inside the liquid stearic acid.

The thermocouple measurements confirm the constant solidification front velocity found previously and give indications about the temperature profile. The temperature profile is flattening with time because the width of the temperature front, w, is increasing, but also because the temperature difference, ∆T = T 1 -T 0 , between the melt and the coolant liquid is slightly decreasing. In fact, the emplacement temperature, T 1 , is decreasing in the experiment from its initial value to the solidification temperature, T S , and consequently

∆T is varying over 15%. This variation is rather small, and does not significantly affect the temperature profile predicted by the conductive cooling model, Eq. ( 3), but seems to be enough to affect the front velocity. While the flattening of the temperature profile should induce a front velocity decrease, the fact that T 1 approaches T S as the experiment proceeds seems to counteract the former effect and maintain a constant solidification front velocity. The value of the temperature front width, w, is determined as the width over which the temperature changes from the coolant temperature, T 0 , to the solidification We interpret the termination of the columnar joints after 2 cm to be caused by the flattening of the temperature profile, such that the temperature gradient falls below a threshold needed for fracturing to occur. We noted the presence of columns with other coolant liquids. Using water at 0 • C, we could observe a columnar structure with about the same column diameter but the columns extended only a few millimeters into the sample.

For the water-cooled samples, the temperature difference, T 1 -T 0 , is about one third of the temperature difference obtained with nitrogen as coolant. The temperature gradient therefore decreases below the threshold value necessary for fracture to occur faster than in the nitrogen cooled case and the columns do not extend as far into the sample.

We noticed that purer samples of stearic acid led to the formation of crystalline structures spreading through the whole sample while it solidified and cooled down. In this case, no columnar structure was observed. Addition of palmitic acid (5-50% in mass), as found in the technical grade, leads to the formation of columns but with a smaller diameter, a better cohesion and with a shiny aspect suggesting the presence of larger crystals. The presence of impurities in the technical grade might prevent the formation of larger crystal structures.

Tensile loading measurements were performed on cylindrical samples (5 mm in diameter, 5 cm in height) prepared as follows: after the wax is melted, we let the liquid stearic acid cool down and solidify slowly at room temperature to prevent the formation of columnar structures. At room temperature, the tensile strength is σ c = (0.34 ± 0.06) MPa and the

Young modulus E = (410 ± 50) MPa.

D R A F T February 6, 2022, 10:10am D R A F T
If a better control of the temperature evolution is gained, we suggest that the stearic acid experiments can be used as a model system for igneous columnar jointing in order to gain insight into the effect of an initial surface crack pattern, to study if two crack fronts approaching each other can produce an entablature like region and to study the effect of the confining geometry. Also the possibility of a range of possible column diameters for each set of material properties and cooling parameters could be investigated. Since we do not yet have the required control of the temperature evolution, the models developed in this paper will instead mainly be compared to numerical simulations and available field measurements for basalt.

A two-dimensional model of columnar jointing

To model the idea of a range of possible column diameters, we consider the idealized version of a system undergoing columnar jointing shown in Figure 3a. We will further simplify the model to an array of semi-infinite cracks in a two-dimensional linear-elastic strip as displayed in Figure 3b. The two-dimensional strip can be thought of as resulting from a plane cut along the red dashed line in Figure 3a. Due to the symmetry of the system, the same plane cut will be encountered repeatedly with a distance of approximately 2 as one moves in the perpendicular Z-direction. We thus expect the Z-direction to be of less importance than the X and Y -directions. Real columnar joint formations are neither two-dimensional, nor perfectly hexagonal, and the validity of the two-dimensional model thus ultimately relies on the agreement between the three-dimensional numerical simulations presented in Section 4 and the predictions of the two-dimensional model.

The two-dimensional strip is displayed in Figure 3b together with a temperature profile, T w,a (X), of the Eq. ( 5) type, which is expected to lead to columnar jointing with a constant D R A F T February 6, 2022, 10:10am D R A F T diameter, . This temperature profile has a fixed shape in a frame of reference moving with the convection zone at speed v as discussed in Section 1.1.2. In contrast to previous sections, we now take the coordinate system (X, Y, Z) to be the frame of reference moving with the convection zone and the crack tips, with the crack tips located at Z = 0. For the purpose of calculating thermal contraction, we are only interested in the temperature deviation from the un-contracted molten interior at infinity and we therefore consider:

T w,a (X) = ∆T 1 -e -(X+a)/w θ(X + a) -∆T (6)
where ∆T = T 1 -T 0 is the maximal temperature difference with respect to the temperature of the undeformed system, a is the signed distance between the temperature front and the crack tips and the length scale, w, describes the width of the temperature profile. T 1 is taken to be the emplacement temperature of the lava and T 0 = 100 • C is the temperature of the steam filled cracks. For columnar jointing in igneous rock, only a > 0 is relevant, since the case a < 0 would imply that the steam/water convection zone penetrates further into the rock than the cracks. We will therefore refer to a as the crack lead length.

In dimensionless form, the temperature field in Eq. ( 6) in the moving frame of reference become:

T P e,δ (x) = T w,a (X) ∆T = 1 -e -(x+δ)P e/2 θ(x + δ) -1 (7)
where we have defined the dimensionless moving frame coordinates:

x = X b y = Y b (8)
and the dimensionless control parameters:

δ = a b P e = 2b w = w = v D (9) D R A F T February 6, 2022, 10:10am D R A F T
where δ is the dimensionless crack lead length and the Péclet number, P e, is the inverse dimensionless width of the temperature profile. For later convenience, we have scaled the parameters with half the column diameter, b, and not the full column diameter, = 2b, as in Eq. ( 5). The only difference generated is a factor of one half in the exponent of Eq. ( 7).

In order to determine the crack growth, we are interested in calculating the mode I stress intensity factor of the crack tips. We therefore need to know how its value depends on the distance between the crack tips and the temperature front δ. We further need to know how it depends on the parameters of the temperature front ∆T and P e as well as on parameters of the solidifying material itself E, α T , and K I,c . Here E is the Young's modulus, α T the coefficient of linear thermal expansion and K I,c the critical mode I stress intensity factor (fracture toughness).

We assume that the cracks propagate straight, such that the mode II stress intensity factor is always zero. Due to the periodicity in the y-direction, we only need to consider the region Y ∈ [-b, b] or equivalently y ∈ [-1, 1], which is colored gray in Figure 3b.

It is also assumed, that the changes in the thermal fluxes happen on a time scale much larger than the time needed for the material to reach elastostatic equilibrium. This is a reasonable assumption since the speed of sound in basalt is of the order 10 3 m/s [START_REF] Gleason | Elasticity of Materials at High Pressure[END_REF] whereas the speed of the thermal front is of the order 10 -8 m/s [START_REF] Hardee | Solidification in Kilauea Iki lava lake[END_REF]. We will not consider the effect of temperature dependent elastic parameters. We define the dimensionless stress, σ ij , and strain, ij : where i, j runs over the dimensionless parameters x, y and I, J refers to the corresponding parameters X, Y with dimension of length. Under plane stress conditions, the elastostatic mechanical problem is described by three sets of equations:

σ ij = σ IJ Eα T ∆T ij = IJ α T ∆T (
Force balance:

∂σ ij ∂x j = 0 (11) xx = (σ xx -νσ yy ) + T P e,δ (x)
Constitutive relations: yy = (σ yy -νσ xx ) + T P e,δ (x)

xy = (1 + ν)σ xy (12)
Compatibility:

2 ∂ 2 xy ∂x∂y = ∂ 2 xx ∂y 2 + ∂ 2 yy ∂x 2 (13)
In order to calculate the stress tensor components, we make use of the dimensionless Airy stress function, φ(x, y), which is related to the stress components through the following derivatives:

σ xx = ∂ 2 φ ∂y 2 , σ yy = ∂ 2 φ ∂x 2 , σ xy = - ∂ 2 φ ∂x∂y . (14)
By inserting Eq. ( 12) in Eq. ( 13) and simplifying using Eq. ( 11), we finally obtain from Eq. ( 14) the equation:

∇ 2 ∇ 2 φ(x, y) = -∇ 2 T P e,δ (x) (15)
The boundary conditions follow from the periodicity of the system and the fact, that σ yy must vanish on the crack faces, whereas the dimensionless displacement in the y-direction, σ yy (x, 1) = 0 for x < 0, u y (x, 1) = 0 for x > 0 (19)

D

Results

The calculation of the stress field from Eq. ( 15) on the entire x-axis is made complicated by the mixed boundary condition in Eq. ( 19) and has to our knowledge not previously been performed. We have used the Wiener-Hopf method along the lines of [START_REF] Marder | Instability of a crack in a heated strip[END_REF] to implement the mixed boundary condition, and solve for the stress intensity factor at

κ I,c = K I,c Eα T ∆T √ b = b min b (21)
then the fracture criteria K I = K I,c can be restated as:

κ I (P e, δ) = κ I,c (22)
and a minimum column spacing, b = b min , is obtained in the case of a uniform contraction

where P e = 0 and the stress intensity factor reaches its maximal value κ I (0, δ) = 1 for all crack lead lengths, δ.

For a given material and temperature profile (fixed b min and w), the fracture criterion in Eq. ( 22 We expect the approximate one-to-one correspondence between P e and κ I,c to be caused by a bifurcation instability [START_REF] Hofmann | Spacing of crack patterns driven by steady-state cooling or drying and influenced by a solidification boundary[END_REF][START_REF] Bahr | Diameter of basalt columns derived from fracture mechanics bifurcation analysis[END_REF]. If we consider all pairs of (P e, δ) fullfilling Eq. ( 22) there will be a minimum value of the Péclet number below which the system is unstable to a small advance of one crack tip relative to its neighbors, leading to a doubling of the column width. It can be motivated, that steady columnar joint formation with a constant width occurs close to the stability limit, by considering the cooling history of the system. When cooling is initiated, the cooling effect of the cracks in the material surface is negligible compared to the external diffusive cooling from e.g. the underlying ground. As the diffusive cooling front slows down with time, a series of bifurcation events are induced as the stability limit for the column width is encountered, and the system coarsens. Eventually, convectional cooling becomes the dominant cooling mechanism, and the temperature field and crack propagation reaches a steady state. The column width will thus approximately correspond to the stability limit.

To obtain a quantitative relationship between the column width and the dimensionless fracture toughness, we will in the following assume a one-to-one correspondence between the Péclet number, P e, and the fracture toughness, κ I,c . This assumption together with Eq. ( 22) implies, that δ, must also be a function of the Péclet number. In addition to the three-dimensional simulations, we have performed simulations in two dimensions, to test more directly the analytical calculations and the relation between P e and δ. In two dimensions, we have measured the average separation between cracks, , and the crack lead length, δ . Figure 6 in agreement with the 2D simulation measurements, see Figure 6, and furthermore yield a best fit between the pairs of κ I,c and measured P e from the three-dimensional simulations, see Figure 7. We note that the data in Figure 6 could be consistent with a narrow but finite range of crack lead lengths, δ, possible for each Péclet number.

Here we assume the power law relation in Eq. ( 24). The power law relation allows us to replace δ with a function of the Péclet number δ = g(P e), and thus analytically obtain a one-to-one relation between the system parameters in κ I,c , and the Péclet number. The fracture criteria in Eq. ( 22) can be restated as:

f (P e) = κ I (P e, g(P e)) = κ I,c (25)
where we have defined the scaling function f (P e) = κ I (P e, g(P e)). The scaling function, plotted in Figure 7, is in excellent agreement with three-dimensional discrete element simulations. This supports the validity of a one-to-one relationship between the Péclet number and the dimensionless fracture toughness, κ I,c . The scaling function agrees reasonably well with estimates of κ I,c and P e for the Kilauea Iki lava lake and stearic acid, where the speed, v, of the temperature front has been directly measured. The estimate of P e and κ I,c for the Kilauea Iki lava lake is based on the parameters listed in Table 1 and a column diameter assumed similar to that of the Prehistoric Makaopuhi lava lake. The Kilauea Iki lava lake is to our knowledge the only published field data, where the speed, v, of the thermal front has been measured directly. The stearic acid estimate of P e and κ I,c is based on the parameters listed in Table 2.

We stress, that the column width predicted by Eq. ( 25) is a function of both the material properties (E, α T , D, K I,c ) and the externally imposed cooling conditions (∆T, v). Figure 8 shows the temperature front speed, v, estimated from our model using the parameters in Table 3 and field measurements of for sites in the Columbia River Basalt Group available from the supplementary material of [START_REF] Goehring | Scaling of columnar joints in basalt[END_REF].

The corresponding estimated Péclet numbers are in the range P e ∈ [0.2 -0.3].

The figure also shows the estimates of v for the same sites based on field measurements of the striae heights, s , as presented in [START_REF] Goehring | Nonequilibrium scale selection mechanism for columnar jointing[END_REF]. The model's estimates, based on , systematically predict values of v which are a factor two smaller than the estimates based on s . The two estimates are close, when keeping in mind that the two methods use independent field data, and s , and different models for the system. We propose our model as a way to quickly estimate the temperature front speed, v, when no striae height measurements are available.

Conclusion

We have used stearic acid as a new model material for columnar jointing caused by thermal contraction. If better temperature control is gained in the stearic acid experiments, we expect that our experimental system could further be used to investigate the entablature formation as well as the possibility of a range of column diameters for each set of material properties and cooling parameters.

D R A F T February 6, 2022, 10:10am D R A F T parameter q = α + iτ , where α and τ are real:

σ+ (q, 1) = ∞ -∞ e iqx σ + (x, 1) dx = ∞ 0 e iqx σ + (x, 1) dx = ∞ 0 e iαx e -τ x σ + (x, 1) dx (A10)
Assuming that the behavior of |σ + (x, 1)| is O(e τσx ) as x → ∞, then the integral in Eq.(A10) converges in the complex half plane given by τ > τ σ . I.e., σ+ (q, 1) is analytic in this half plane.

If σ + (x, 1) has algebraic growth, is constant or decays slower than an exponential function for x → ∞, then τ σ = 0. So the real axis is only a part of the region of analyticity if σ + (x, 1) has exponential decay at infinity.

We can define ũ-(q, 1) by doing similar considerations. We know that u -(x, 1) vanishes in the uncracked material for x > 0:

ũ-(q, 1) = ∞ -∞ e iqx u -(x, 1) dx = 0 -∞ e iαx e -τ x u -(x, 1) dx (A11)
Assuming the behavior of |u -(x, 1)| is O(e τux ) as x → -∞, then the integral in Eq. (A11) converges in the complex half plane given by τ < τ u , such that ũ-(q, 1) is analytic in this region.

Solving Eq. ( A9) with the boundary conditions of Eq. (A2-A4, A8) yields the Wiener-Hopf Eq. (A12):

σ+ (q, 1) = D P e,δ (q) + F (q)ũ -(q, 1)

(A12) D R A F T February 6, 2022, 10:10am D R A F T
where ũ-(q, 1) is still partly unknown. The functions F (q) and D P e,δ (q) are defined as:

F (q) = q 2 coth(q) + q csch 2 (q) (A13)
D P e,δ (q) = -TPe,δ (q) (A14)

A2. Solving the Wiener-Hopf Equation

The Wiener-Hopf method allows us, to use our knowledge of the asymptotic behavior of u -(x, 1) and σ + (x, 1) to express σ+ (q, 1) in Eq. ( A12) purely in terms of the known functions F (q) and D P e,δ (q) and without reference to the partly unknown u -(x, 1).

The first step is to separate Eq. ( A12) into a right hand side analytic in a half plane τ < τ u and a left hand side analytic in a half plane τ > τ σ . Provided that the two half planes overlap, i.e. τ σ < τ u , then the right hand side represents the analytic continuation of the left hand side into the entire complex plane and vice versa.

To accomplish the separation of Eq. ( A12), we need to factorize the function F (q) into two parts that are analytic in each their respective half plane of the complex plane. The function F (q) has poles on the imaginary axis for q = iπn with n ∈ Z/{0} and increases linearly as q/2 for |q| → ∞.

We choose to factorize F (q) into a fraction F (q) = F -(q)/F + (q) where F -(q) has no poles or zeros for τ < τ F -and F + (q) has no poles or zeros for τ > τ F + , where we again use the complex Fourier variable q = α + iτ . The factorization was performed numerically using 32 Chebyshev polynomials following the strategy of [START_REF] Liu | The energy of a steady-state crack in a strip[END_REF]. The large q behavior of F + and F -was chosen to be:

lim |q|→∞ F + (q) ∼ 2 -iq (A15) lim |q|→∞ F -(q) ∼ iq 2 (A16) D R A F T February 6, 2022, 10:10am D R A F T
which produces the correct behavior lim |q|→∞ F (q) ∼ q/2 for large |q|. Using the factorized F (q) we can rewrite the Wiener-Hopf Eq. (A12) as:

σ+ (q, 1)F + (q) = D P e,δ (q)F + (q) + F -(q)ũ -(q, 1) = e -iqδ D P e,0 (q)F + (q) + F -(q)ũ -(q, 1) (A17)

where we have used the translation properties of the Fourier transform to pull the exponential e -iqδ out in the last line. We need to separate the product D P e,0 (q)F + (q) into regions of analyticity. We use the Fourier convolution theorem:

D P e,0 (q)F + (q) = ∞ -∞ e iqx dx ∞ -∞ D P e,0 (x )F + (x -x )dx = ∞ -∞ e iqx h P e (x) dx = 0 -∞ e iqx h P e (x) dx + ∞ 0 e iqx h P e (x) dx = h -(q) + h + (q) (A18)
Where h + (q) does not have any poles in the upper τ > τ h+ half of the complex plane and h -(q) does not have any poles in the lower τ < τ h-half of the complex plane. With this factorization, the Wiener-Hopf Eq. (A17) takes the form:

σ+ (q, 1)F + (q) -h + (q)e -iqδ = h -(q)e -iqδ + F -(q)ũ -(q, 1) (A19)
The left hand side is analytic in the upper half plane τ > max(τ σ , τ F + , τ f + ), and the right hand side is analytic in a lower half plane τ < min(τ u , τ F -, τ f -). Provided that the right

σ+ (q, 1)F + (q) -h + (q)e -iqδ = h -(q)e -iqδ + F -(q)ũ -(q, 1) = W (q) (A20)
We expect that u y (x, b) goes to a constant for x → ∞ which means that τ u = 0. There will be an overlap of the two regions of analyticity if σ yy (x, 1) has just a slight exponential decay for x → ∞, such that τ σ = -where is a small positive real number. The overlap is depicted as the gray shaded area in Figure 9 and an overview of the regions of analyticity is shown in Table 4. It is not physically unreasonable to expect a slight exponential decay of σ yy (x, 1) as x → ∞, since we know σ yy (x, 1) does decay in this limit.

Assuming τ σ = -, the function W (q) is entire. To determine the constant value of W (q) we use Abel theorems, which relate the asymptotic behaviors of a function g(r) and its one sided Fourier transforms g+ (q) or g-(q). We will use the following Abel theorem [Freund , 1990, Sec. 1.3] rewritten for a plus type Fourier transform with q = α +iτ , where both α and τ are real:

lim x→0 + Γ(1 + γ)x -γ g(x) = lim τ →+∞ (-iq) 1+γ g + (q) (A21)
Which is valid for γ > -1. We can extract the asymptotic behavior of the Fourier transformed stress function using the Abel theorem:

lim x→0 + σ yy (x, 1) ∼ lim x→0 + 1 √ x ⇒ lim τ →+∞ σ + (q, 1) ∼ lim τ →+∞ 1 √ -iq (A22)
From the behavior of F + (q) and D P e,0 (q) we directly have:

lim |q|→∞ h + (q) ∼ lim |q|→∞ 2 -iq 1 iq (A23) D R A F T February 6, 2022, 10:10am D R A F T
If we consider the left hand side of Eq. (A20), using the asymptotic results of Eq. (A15, A22, A23) we get:

lim τ →+∞ W (q) = lim τ →+∞ (σ + (q, 1)F + (q) -h + (q)) ∼ lim τ →+∞ 1 √ -iq • - 2 iq - 2 -iq 1 iq = lim τ →+∞ √ 2 τ -iα + 2 τ -iα 1 τ -iα = 0 (A24)
Now, W (q) is analytic in the entire plane, bounded and equal to zero at infinity. Therefore, by Liouvilles Theorem, it must be identical zero in the entire plane. This allow us to finally determine σ+ (q, 1) from Eq. (A20) without any reference to the partly unknown ũ-(q, 1):

σ+ (q, 1) = 1 F + (q) h + (q)e -iqδ = 1 F + (q) ∞ 0 dx e iqx h P e (x + δ) = 1 F + (q) ∞ -∞ dx e iqx θ(x)h P e (x + δ) (A25)
We note, that from Eq. (A18) we can also express h P e (x) as:

h P e (x) = ∞ -∞ dq 2π e -iqx F + (q)D P e,0 (q) (A26)
where the inverse Fourier transform integral has to be performed within the strip of analyticity of the integrand.

A3. Extracting the stress intensity factor

The Abel theorems in Eq. (A21) allow us to infer the x → 0 + behavior of σ yy (x, 1) from the τ → ∞ behavior of σ+ (q, 1), where again q = α + iτ .

D R A F T February 6, 2022, 10:10am D R A F T
We define the dimensionless mode I stress intensity factor κ I as:

κ I (δ) = lim x→0 + √ 2π x 1/2 σ yy (x, 1) (A27) = lim τ →+∞ √ 2π Γ 1 2 (-iq) 1/2 σ+ (q, 1) (A28) = lim τ →+∞ -2iq σ+ (q, 1) (A29)
To go from Eq. (A27) to Eq. (A28), we have used the Abel theorem in Eq. (A21) for γ = -1/2. From Eq. (A25), we see that the large q behavior of σ+ (q, 1) is partly due to the 1/F + (q) factor and partly due to the integrand. The large q behavior of the integrand in Eq. ( A25) is governed by the step function discontinuity of the integrand at x = 0.

The large q behavior of the integrand becomes:

lim |q|→∞ ∞ -∞ dx e iqx θ(x)h P e (x + δ) = lim |q|→∞ -1 iq h P e (δ) (A30)
where the Fourier transform of the step function θ(x) is -1/iq and exists for τ > 0. The large q behavior of F + (q) we already know from Eq. (A15). Collecting our knowledge on σ+ (q, 1) we arrive at:

κ I (δ) = lim τ →+∞ -2iq -iq 2 -1 iq h P e (δ) = lim τ →+∞ 2(τ -iα) τ -iα 2 1 τ -iα h P e (δ) = h P e (δ) (A31)
This is the result stated in Eq. ( 20) of the main text.

with their nearest-neighbors (NN) and next-nearest-neighbors (NNN) through Hookean springs with spring constants k 1 and k 2 for NNs and NNNs, respectively. The spring connecting the nodes i and j has an equilibrium length L eq ij and holds an elastic energy quadratic in the deviation from the equilibrium length,

U ij = 1 2 k ij (L ij -L eq ij ) 2 (B1)
where x i denotes the position of the i'th node and L ij = |x i -x j |. The total energy of the system is given by the sum over all springs U tot = U ij . The elasticity tensor corresponding to these interactions can be written in the form [START_REF] Hayakawa | Pattern selection of multicrack propagation in quenched crystals[END_REF]:

C iiii = (k 1 + 2k 2 )/L 0 (B2) C iijj = C ijji = C ijij = k 2 /L 0 . (B3)
The thermal contraction of the system is modeled by changing the local equilibrium distance between the nodes. Like in the previous section, we shall consider a contraction front on the form of Eq. ( 6):

L eq ij (X, X 0 ) = L 0 ij 1 + α T ∆T 1 -e -X-X 0 w θ (X -X 0 ) -α T ∆T (B4)
where X 0 = vt = D w t is the position of the temperature front, which moves through the system with time. Since elastic relaxation happens many orders of magnitude faster than thermal relaxation, the simulations proceed by advancing the contraction front a length dX 0 , relaxing the system to elastostatic equilibrium, and repeating the process.

The value of dX 0 does not affect the resulting fracture pattern. The state of elastostatic equilibrium or force balance is reached by minimizing the elastic energy, i.e. finding the lattice configuration for which ∇U = ∇ U ij = 0. To that end, we apply a conjugate D R A F T February 6, 2022, 10:10am D R A F T gradient method in the multidimensional configuration space. In the simulations, we use a failure criterion where an elastic bond in the lattice is broken when it is strained, ij , beyond a critical strain c :

ij > c . (B5)
With this threshold, bonds do not break during compression, i.e. the fractures forming the columnar joints will mostly grow at the trailing edge of the contraction front. Note, that while this asymmetry in the fracture criterion have little impact on the morphology of the steady state patterns, it may influence the fracture patterns formed close to the boundary of the system, i.e. where the front is initiated.

If more than one bond is strained beyond the critical level, the bond with the highest strain is broken and the lattice is relaxed to the new state of elastostatic equilibrium. We then check whether other bonds still exceed the critical stress. If a bond still exceeds the critical stress the process is repeated until no such bond exists and then finally we advance the contraction front by a small increment, dX 0 . In that way, bonds are broken one-by-one until none of them exceeds the critical strain.

Note that in the simulations the temperature profile starts out with the steady state shape and therefore the initial cracks nucleating at the surface of the simulation system might not have any resemblance to the columns close to the cooling surface in real system. it from the computational mesh (by labeling it as air). This is followed by reassembling the stiffness matrix and re-equilibrating the system. Those two steps are then repeated as long as there are elements with σ M > σ c in the mesh. Re-equilibration, which amounts to solving the constitutive equations with a new stiffness matrix, is a computationally intensive process hence, in order to improve the method's performance we introduce two simplifications:

• System relaxation in the fracture resolution step is performed on a submesh composed of elements lying in the X-distance 0.1 from the temperature front (which is equivalent to treating the excluded region as fluid).

• Instead of only removing a single element at a time, all elements whose von Mises stress are greater than γσ c , γ > 1 are all removed at once, where γ is chosen in a way that it does not affect the resulting fracture pattern (determined in our simulation to be γ = 1.1).

The resulting set of removed tetrahedra displayed in Figure 10 The values are all taken from 3: this paper, section 2. The temperature front width w is from thermocouple measurements.

b No measurements of the Poisson's ratio of stearic acid were found in the literature, The common strip of analyticity is gray shaded and covers the region -< τ < 0 in the complex plane q = α + iτ .

Table 4. Overview of the regions where the +/-functions are analytical and nonzero. The zeros of F (q) closest to the real axis were numerically determined to be at τ = ±2.1062. Thus, the zeros occur closer to the real line than the first poles at τ = ±π, and the value of τ F -and τ F + are set by these zeros.

Upper half plane Lower half plane

τ σ - τ u 0 τ h + -2.1062 τ h - 0 τ F + -2.1062 τ F - 2

  temperature, T S . The value of w increases with time and we estimate its range of values as w = (4 ± 2) mm for our thermocouple measurements. The experimental values for D R A F T February 6, 2022, 10:10am D R A F T the average column diameter and front width are used to estimate a Péclet number of P e = /w = 0.5 ± 0.4 for the stearic acid system.

  10) D R A F T February 6, 2022, 10:10am D R A F T

  R A F T February 6, 2022, 10:10am D R A F T u y = u Y /b, must vanish on the rest of the crack line. σ xy (x, 1) = σ xy (x, -1)

  the crack tip. The solution is derived in Appendix A and we obtain the dimensionless mode I stress intensity factor κ I = K I /(Eα T ∆T √ b): κ I (P e, δ) = h P e (δ) (20) where δ = a/b is the dimensionless crack lead length and the Péclet number, P e = 2b/w, is the inverse dimensionless width of the temperature profile. The function h P e (δ) is found analytically in the Appendix, Eq. (A26), and the stress intensity factor is plotted for different values of P e and δ in Figure 4. Crack advance in columnar jointing basalt occurs incrementally, producing striae as mentioned in the introduction. However, we will not try to model this feature, but consider the crack propagation to be in steady state, such that K I = K I,c always holds at the crack tip, where K I,c is the material dependent critical mode I stress intensity factor (fracture toughness). The crack thus propagates continuously with the same speed as the D R A F T February 6, 2022, 10:10am D R A F T temperature front. If we define the dimensionless critical stress intensity factor, κ I,c , and the mechanical loading length, b min , as:

  ) can be fulfilled for all b ≥ b min by changing the crack tip position, δ = a/b, and thus adjusting κ I (P e, δ) in the interval [0; 1]. The physical mechanism behind the continuous set of possible (P e, δ)-pairs for fixed material and temperature parameters, can be thought of as follows: If the fracture spacing, and thus the Peclét, P e, is imagined to increase slightly, then the resulting fracture density decreases and the stress intensity at the crack tips increases. To accommodate this change, the crack tips can move further in to the hot un-contracted region by increasing δ and thus bring the stress intensity at the crack tips back to the original level.In summary, the two-dimensional model of an infinite array of cracks predicts, that there is a minimum column diameter for a given material and temperature profile: b ≥ b min (23) where b min = (K I,c /Eα T ∆T ) 2 and the range of b-values fulfilling Eq. (23) are equally possible within the framework of the model. D R A F T February 6, 2022, 10:10am D R A F TTo test whether a range of column diameters occurs, we have performed a set of discrete element simulations (Appendix B). In the simulations the temperature front in Eq. (6) is propagated through a three dimensional linear elastic material. The material is represented by discrete elements, i.e. connected springs which thermally can contract and which can break if a critical stress or strain is reached. An example of the resulting fracture network is shown in Figure5. The dimensionless fracture toughness, κ I,c , and the temperature front width, w, are simulation parameters, whereas the column diameter, , is measured when the system reaches a state of steady crack propagation with constant column diameter. From the measured column diameter, the Péclet number of the simulation, P e = /w, can be found.We have performed simulations where the width of the temperature front, w, is either suddenly or slowly increased/decreased for fixed κ I,c . In general, we observe a corresponding increase/decrease in the column diameter such that the Péclet number stays approximately constant. If all column diameters in the range b ≥ b min were equally possible, we would expect to see instances where the column width remains constant at the cost of a change in distance between the temperature profile and the crack tips. As w changes, but b stays fixed, the result is a change of the Péclet number. In summary, the simulated systems seem to select one specific column diameter, , for each temperature profile width, w, leading to a one-to-one correspondence between the Péclet number, P e, and the dimensionless fracture toughness, κ I,c . However, the possibility of a narrow but finite range of allowed Péclet numbers, and thus column widths, can not be ruled out from the present simulations. D R A F T February 6, 2022, 10:10am D R A F T

  indicates a power law relation between the Péclet number, P e, and the crack lead length, δ, measured in two dimensions: δ = g(P e) = c 1 P e c 2 (24) D R A F T February 6, 2022, 10:10am D R A F T As the error-bars are large, we propose the coefficients c 1 = 5.4 and c 2 = -1.4 which are

D

  R A F T February 6, 2022, 10:10am D R A F T 4.1. Prediction of the temperature front propagation speed from the scaling function The scaling function, Eq. (25), allows us to estimate the Péclet number of a process P e = f -1 (κ I,c ), from the dimensionless fracture toughness κ I,c = K I,c /Eα T ∆T √ b. From the Péclet number, P e = /w = v /D, the temperature front propagation speed, v, can be estimated, if the thermal diffusivity, D, and the column diameter, = 2b, is known.

  hand side and the left hand side overlap in some non-empty open subset like a line or a small area, they are each others analytical continuation into the other half plane and D R A F T February 6, 2022, 10:10am D R A F T equal to the same entire function W (q):

  After a transient dynamics over a length of approximately 2w, the cracks in the simulations form a regular polygonal pattern equivalent to the columnar joints. An example of a simulation is shown in Figure 5. When the regular polygonal crack pattern has formed, then the average column diameter is measured. From this measurement, the Péclet number P e = /w for the simulation can be determined. D R A F T February 6, 2022, 10:10am D R A F T time step advances the contraction front a step dX 0 , updates the equilibrium von Mises stress values, and proceeds to fracture resolution. The fracture resolution is performed by finding the element with the maximum von Mises stress exceeding σ c and removing

Figure 2 .Figure 4 .Figure 5 .Figure 6 .Figure 7 .

 24567 Figure 1. a) Exposed columnar joint formation at Svartifoss, Iceland. b) Sketch

a

  The values are taken from 1:[START_REF] Degraff | Effect of thermal regime on growth increment and spacing of contraction joints in basaltic lava[END_REF], 2: Ryan and Sammis[1981] and 3:[START_REF] Hardee | Solidification in Kilauea Iki lava lake[END_REF].b A 10% deviation of E is inferred from Figure2of[START_REF] Ryan | The glass transition in basalt[END_REF].cis the column diameter of the Prehistoric Makaopuhi lava lake, which is assumed similar to the Kilauea Iki lava lake. The average column face width has been convertedto diameter using a hexagonal column shape.D R A F T February 6, 2022, 10:10am D R A F T

Figure 8 .

 8 Figure 8. The temperature front propagation speed, v = P e D/ , is estimated for dif-

[Figure 9 .

 9 Figure 9. Regions of analyticity in the complex plane for different functions. The left

Figure 10 .

 10 Figure 10. Example of the resulting fracture network in a finite element simulation of

Table 1 .

 1 Physical properties of basalt and parameters of the temperature front at the Kilauea Iki lava lake a

	Property	Symbol Value		Source
	Fracture toughness	K I,c 2.2 MPa √	m	1
	Young's modulus	E	(54.5 ± 5.5) GPa	2 b

Table 2 .

 2 Physical properties of stearic acid and parameters of the temperature front a

	Property	Symbol Value	Source
	Tensile strength	σ c	(0.34 ± 0.06) MPa	3
	Young's modulus	E	(410 ± 50) MPa	3
	Poisson's ratio	ν	0.4	3 b
	Linear contraction	α T ∆T 1 3 10%	3
	Temperature front width	w	(4 ± 2) mm	3
	Column diameter		(2 ± 1) mm	3
	a

Table 3 .

 3 Physical properties of basalt and parameters of the temperature front in the Columbia River Basalt Group a

	Property	Symbol Value	Source
	Fracture toughness	K I,c 2.2 MPa √ m	1
	Young's modulus	E	50 GPa	5
	Poisson's ratio	ν	0.2	5
	Linear coefficient	α T	7 • 10 -6 / • C	5
	of thermal expansion
	Thermal diffusivity	D	6.5 • 10 -7 m 2 /s	5
	Max. temperature difference ∆T	990 • C	5

a The values are taken from 1: Degraff and Aydin

[1993]

and 5: Goehring and Morris

Acknowledgments. This study was supported by the grant Earth Patterns by Villum Fonden. Fabrice Mortessagne is thanked for his support of the project and the CNRS is thanked for financial support. The data for this paper are available by contacting the corresponding author at mathies@nbi.dk.

Based on numerical simulations and analytical calculations, we argued that one Péclet number is selected for each dimensionless critical stress intensity factor, though a narrow finite range of Péclet numbers could not be ruled out.

We derived the functional form of the relation between the Péclet number and the dimensionless critical stress intensity factor and finally showed that our derivations are in reasonable agreement with data obtained from stearic acid experiments, three-dimensional numerical simulations and field data. We suggest our scaling function as a quick way to determine the rate by which the fracture front and therefore the cooling front advances in a system using as input basic properties of the material, the emplacement temperature and field measurements of the column diameter. Alternatively, an estimate of the column diameter can be obtained from the scaling function if the fracture front speed is known.

Appendix A: Analytical solution for the stress intensity factor

A1. The Wiener-Hopf Equation

We want to solve Eq. ( 15), which is defined on the entire x-axis with the boundary conditions in Eq. (16)(17)(18)(19). For convenience, Eq. ( 15) is here reproduced as Eq. (A1):

and the boundary conditions are reproduced as Eq. (A2-A5):

where x, y are dimensionless coordinates, σ ij , ij are dimensionless stress and strain,

T P e,δ (x) defined in Eq. ( 7) is the dimensionless temperature front, δ is the dimensionless distance between the crack tips and the temperature front and P e is the dimensionless width of the temperature front.

We note, that had we instead of plane stress worked with plain strain conditions, the modification of the constitutive relations would result in a factor of 1/(1 -ν 2 ) multiplying the right hand side of Eq. (A1). All calculations to be performed in this section can therefore easily be modified to a plane strain setting.

If we introduce the partly unknown functions u -(x, 1) and σ + (x, 1):

then the boundary condition for u y (x, 1) in Eq. ( A5) can be extended to the entire x-axis and restated as:

where u -(x, 1) is partly unknown so far.

Eq. (A1) can be Fourier transformed in the x-direction and turned into a fourth order inhomogeneous ODE. We denote the dimensionless transform variable q and the transformed functions by a tilde:

We now consider the Fourier transform of the functions in Eq. (A6-A7) utilizing that they both vanish on half of the x-axis. For the Fourier transform of σ + (x, 1) we use a complex applied on the faces parallel to the X-axis. The computational domain is finely discretized using tetrahedral elements (4 million tetrahedra). We assume the material to be isotropic and Cauchy elastic in the simulation, with the Young's modulus E and Poisson ratio ν = 0.25. The assumption of small displacements is valid in the context of columnar joint formation and the elasticity model is appropriate for simulating infinite systems.

We use the von Mises yield criterion, which states that failure occurs when the von Mises stress, σ M , exceeds the material's yield strength, σ c , chosen in simulation to be equal σ c = 0.01E. We use a standard, linear Galerkin method to discretize the constitutive equation.

The propagating temperature front and corresponding thermal contraction are modeled as a body force:

where the maximum contraction is α T ∆T = 0.2, the contraction front is located at X 0 = vt = D w t and the slope 1/w varies between 12 and 18. The simulation method at each D R A F T February 6, 2022, 10:10am

Table 5. List of simulation parameters used in the discrete element simulations.

Simulation parameter Value

Poisson's ratio ν = 0.2 a In Figure 7 the value of κ I,c has been multiplied with (1 -ν) to make it comparable with the scaling function derived for 2D.

Critical strain

Table 6. List of simulation parameters used in the finite element simulations.

Simulation parameter Value

Poisson's ratio ν = 0.25