Geometric Attributes of Polygonal Crack Patterns in Columnar Joints

Abstract

A well‐ordered polygonal crack pattern is frequently observed on the outcrop surface of columnar joints. The crack pattern's geometry and topology result from the volume contraction of solidified lava during thermal dissipation. In this study, we conducted field measurements of columnar joints at four different sites with various lithologies in Japan to understand the influence of geography and lithological character on their geometric characteristics. Drone‐aided photography and computer‐based image analysis enabled us to examine the constituent polygons' geometric attributes and their statistics with high accuracy. This methodology revealed a relationship between the geometry and topology of the polygonal crack patterns. Moreover, we established a positive correlation between the polygon area and the number of vertices of a polygon at each study site, irrespective of lithology, or geographic position. These results are consistent with previous investigational studies, in particular starch‐based analog experiments.

Introduction

Columnar jointing is a spectacular geological structure composed of many polygonal prisms of igneous rock or ignimbrite. Formation of columnar jointing is considered a spontaneous process in which solidified lava's rapid cooling causes volume contraction, followed by cracking at the surface. Through further volume contraction, surface cracks propagate toward the inside of the rock, resulting in an array of columns with polygonal cross-sections (Mallet, 1875).

These regularly ordered columnar structures have attracted attention for more than 300 years. Sir R. Bulkeley published the first report on columnar jointing in 1693 on the Giant's Causeway in Northern Ireland, mentioning the frequency with which hexagonal and pentagonal columns of ∼18–20 inches in diameter are found (Bulkeley, 1693). Similar trends in the frequency of polygonal shapes have been observed in other columnar joint localities around the world, including the Devils Postpile in California and Devils Tower in Wyoming (Beard, 1959; Budkewitsch & Robin, 1994; Hetényi et al., 2012). In addition, the observed column diameter is known to vary depending on the location. Fieldwork and crystal texture analysis (Long & Wood, 1986) have suggested that the original lava's cooling rate determines the typical column diameter at a specified location, whereby slower cooling rates produce large-diameter columns and vice-versa. This theory was reconfirmed by a starch-based analog experiment, which proved that the intercrack spacing is inversely proportional to the cracking advance speed (Goehring et al., 2009; Toramaru & Matsumoto, 2004). These findings account for typical column diameters' location-dependence, implying that the columnar diameter is uniform to a considerable degree within a specified location.

From a theoretical viewpoint, a plausible origin of the regularity in the typical column diameter within a specified location can be described by the coupled conductive and convective heat transfer model (Budkewitsch & Robin, 1994; DeGraff & Aydin, 1987, 1993), as shown in Figure 1. The model assumes that upon cooling of lava, internal heat is first dissipated by heat conduction. The conduction-driven heat dissipation causes the solidification and volume contraction of lava such that, beyond the threshold of contraction, the surface cooling fractures and crack fronts propagate inward with time (Grossenbacher & McDuffie, 1995). Within the narrow gap of the cracks, water-steam convection promotes heat transfer, both from the interior of columns and below the crack fronts (Hardee, 1980), eliciting a steep thermal gradient near the crack tip. Guided by the thermal gradient, crack fronts cut into the more massive column (Column B in Figure 1), which reduces the difference in the diameters of the two adjacent columns separated by the narrow gap. Repetition of this process, driven by conduction through lava and convection in the opening, yields a regularly ordered columnar structure with limited diameter variation, as partly proved in numerical simulations (Bahr et al., 2009; Goehring & Morris, 2006; Jagla, 2004; Jagla & Rojo, 2002; Nishimoto et al., 2007; Tang et al., 2017). These joint increments propagate inward and can be observed as stria, which resemble chisel-marks (DeGraff & Aydin, 1987, 1993; Goehring & Morris, 2008; Piombo & Dragoni, 2018). Stria form perpendicular to the direction of crack motion in columnar joints and record the stick-slip type of individual cracks growth. This scenario suggests that columnar joints are distinct geometric patterns created by volcanic activity and subsequent thermal dissipation.

Japan is one of the world's most volcanically active countries, with 111 active volcanoes (List of active volcanoes in Japan, 2018). The significant number of volcanoes is a direct result of Japan's location on the boundary of four tectonic plates: the Pacific Plate, the North American Plate, the Eurasian Plate, and the Philippine Sea Plate. Many volcanoes are distributed along the subduction zones of these plates; consequently, various columnar joints, showing diversity in lithological character and columnar cross-section size, are distributed at more than 60 locations in Japan. Most notably, the “Tatami-Ishi” in southernmost Japan's Okinawa prefecture is a large-scale andesite columnar joint complex extending approximately 600 m off the coast. A salient feature of Tatami-Ishi is a flat exposed surface, which resembles an artificial stone pavement, due to repeated erosion by sea waves. Additionally, the polygonal cracking pattern developed in the rock surface can be observed only at low tide, as the body is submerged at high tide. Such columnar joints that spread over a flat, vast space are infrequent worldwide. In addition to Tatami-Ishi, several other columnar joints occur in Japan, and all show well-ordered polygonal cracking patterns at the exposed flat surfaces. Investigating the geometric features of such cracking networks at each location is essential for obtaining insight into the effects of lithological and geographical conditions on columnar joint formation. Despite an academic interest, however, investigations on columnar joints' geometric features in Japan are limited.

This study conducted field surveys of several columnar joint complexes in igneous rocks of various compositions in Japan. At each locality, the columnar jointing's exposed surface was photographed from above by a drone, followed by image analysis to study the polygonal cracking network's statistical properties. In the image analysis, we focused on three geometric attributes: the polygon area, the number of vertices, and the interior angle at vertices (Figure 2); precise definitions of the three features are given below. A possible relationship between the geometric attributes was also examined.

Geological Properties of Columnar Jointing Localities

Our first task was to obtain aerial photographs of the exposed surfaces of columnar joints. It was necessary to choose field survey sites that did not have significant obstructions, such as buildings or overgrown tree branches, to obtain sufficiently clear image analysis photographs. It was also preferred that the exposed surface did not display significant terrestrial undulation in height, such as that exhibited by a cliff because terrestrial undulation may cause shadows and positional deviations in the aerial photographs that hinder the accurate evaluation of the polygon geometry. Based on these criteria, we selected four sites: Tatami-Ishi in Okinawa (O-site), Tatami-ga-Fuchi in Yamaguchi (Y-site), Hino-Misaki in Shimane (Sm-site), and Tawara-Iso in Shizuoka (Sz-site).

Figure 3 indicates the locations of the four sites on a map of Japan. Three areas (O, Sm, and Sz) are positioned along the seacoast, and the remaining site (Y) is situated inland. Table 1 summarizes the GPS coordinates and lithological character of each locality. Note that the Y-site has the fewest polygons (890) among the four localities. Given the statistical significance of our results discussed below, this value indicates that a sufficient number of polygons were present at each site to determine the cracking patterns' statistical properties. Figures 4a–4d represent the geological maps for each location.

Table 1 GPS Coordinates and Lithologic Character of the Study Sites

Location	Latitude	Longitude	Lithologic character
Tatami-Ishi (O-Site)	N26°20′8.79″	N126°49′27.93″	Andesite
Tatami-Ga-Fuchi (Y-Site)	N34°32′52.65″	N131°40′31.06″	Basalt-basaltic andesite (alkaline basalt, trachybasalt)
Hino-Misaki (Sm-Site)	N35°26′4.09″	N132°37′50.24”	Rhyolite
Tawara-Iso (Sz-Site)	N34°39′32.15″	N138°59′3.51″	Andesite

4 Figure

Geological map of study sites indicated by red open circles: (a) O-site, (b) Y-site, (c) Sm-site, and (d) Sz-site. Right figures are reprints from the Seamless Digital Geological Map of Japan (2015). Left figures are reprints from the Terrain with Labels (2016).

The O-site (Figures 4a and 5) is located along the southeast coast of Kumejima Island, one of Japan's southernmost remote islands. The island is composed mainly of andesite rocks. At O-site, columnar jointing occurs in the Late Miocene to Pliocene volcanic rocks. The columnar joints' planar surface extends along the seacoast for 110 m in a north-south orientation and 30 m east-west. We estimated the typical column diameter at approximately 2 m, considerably more extensive than those found at the other sites. As shown in Figure 5a, a noteworthy observation at O-site is the centimeter-high upward swelling of each polygon's edges and center. Among our four localities, this swelling feature was found only at O-site.

At Y-site (Figures 4b and 6), the cracking network extends for over 2 km along the Tama River. The river basin is a caldera structure formed by the ancient volcanic activity of Mt. Irao, part of the Abu monogenetic volcano group, covering 400 km² and consisting of basaltic to dacitic lava flows. The terrain around this area is composed of Pleistocene-aged lava (ca. 126,000 years old) that flowed into and buried the Tama River's valleys and was subsequently cut down by the erosive activity of the river to form a small gorge. Columnar jointing at Y-site is thought to be developed in alkaline basalt or trachybasalt, while a radioactive dating experiment indicated that Y-site consists of basalt-basaltic andesite (Kakubuchi et al., 2000). Here, the typical column diameter is ∼0.3 m, and column heights exceed 8 m.

At Sm-site (Figures 4c and 7), regular arrays of subaerially exposed, minute columns ∼0.05–0.10 m in diameter are observed along the sea cliffs facing the Sea of Japan. These columns, composed of rhyolite, erupted approximately 16 Ma; the rhyolite's flat surface is thought to be a coastal terrace formed in the Pleistocene and extends for over 5 km to the east. The columns' uplift height is ∼100 m, and the terrace surface descends toward the sea to ∼20 m above sea level. It is interesting to note that the top surface of the Sm-site exhibits distinct long and nearly straight cracks, which look similar to master joint (or primary cracks) reported in the earlier work (Peck & Minakami, 1968; Spry, 1962). These cracks were found only at the Sm-site among the investigation sites we examined. In columnar joint formation, primary cracks are known to form prior to the appearance of other short cracks that constitute the polygonal pattern. Then, as the short cracks grow, they reach the existing primary cracks perpendicularly, causing the so-called “T”-junction (with a combination of joint angles: 90°–90°–180°). The presence of many T-junctions at the Sm-site, as clearly shown in Figures 7b and 7c, is therefore considered to be a remanent of the primary cracks formed in the early stage of columnar joint formation.

At Sz-site (Figures 4d and 8), columnar jointing developed within an andesitic sill of Late Miocene to Pliocene age. Uplift and erosion following the collision of a series of volcanic ranges on the Philippine Sea Plate with the island of Honshu resulted in the exposure of columnar joints on the former seabed. The size and shape of columns observed at Sz-site are comparably uniform, with a typical column diameter of ∼0.2 m.

Methods: Aerial Photography and Image Analysis

Aerial images of the localities were taken using a drone (DJI, Phantom4 Pro+). The camera resolution was set to 1,280 × 720 (high definition), which is sufficient to analyze the geometrical properties of the polygonal crack patterns. First, the drone was directed to rise to a height of 10 m. Next, the camera was pointed at the ground, and a bird's-eye view was taken in hovering mode. The drone was then moved slightly parallel to the ground such that additional overlapping photos could be acquired. This “shoot-and-move” step was repeated until the entire crack pattern at a locality was obtained. The total number of images acquired was 752 at O-site, 826 at Y-site, 209 at Sm-site, and 310 at Sz-site.

Photos taken at each site were merged into a single large-scale image using PhotoScan Professional 1.2.6 (Agisoft) software. While processing the image, we applied orthorectification to remove distortion caused by tile or terrain relief. Orthorectification is an image transformation technique through which an aerial photograph's central perspective is altered to an orthogonal view of the ground. The resultant orthorectified image has a constant spatial scale free from distortion, allowing for an accurate evaluation of distances, angles, and areas associated with the observed polygonal crack patterns.

Following orthorectification, we converted the images into polygonal patterns using ArcGIS (Esri), a geometric tool. This conversion procedure is illustrated in Figure 9. First, all cracks displayed in the orthorectified image (Figure 9a) were manually traced using points and line segments (Figure 9b). Points indicate vertices of polygons, whereas line segments delineate the edges. Manual tracing was necessary as it is difficult for the ArcGIS software to identify any blurred cracks, mainly when the colors of cracks and rocks are relatively similar. We then eliminated unclear polygons that, for example, were covered with sand or seawater. The resulting data set of traced points and line segments, called a shapefile, contains geometric information about the point's coordinates and the topology of the lines joining the points. Next, the area enclosed by line segments was extracted from the shapefile (Figure 9c) and is defined by the particular cross-sectional area of columns. Afterward, the number of vertices and the interior angles at vertices were calculated using point coordinates' data. Once the numeric data were processed, we calculated the three geometric attributes' probability distributions to obtain insight into their statistical properties and intercorrelation. To the best of our knowledge, no previous attempts have been made to measure the three geometric quantities at the four sites to reveal their probability distribution curves.

Results

Polygonal crack pattern geometry extracted from the aerial photos was characterized by the following three parameters: the polygon area (i.e., the cross-sectional area of a column), the number of vertices, and the joint angle at a crack intersection. Figure 2b schematically illustrates these attributes. We computed the geometric attributes of all the polygons contained in the large orthorectified images for all four sites. We then derived the probability distribution of the geometric attributes to determine the similarities and differences in the localities' statistical properties.

Statistics of Geometric Attributes

Figures 10a–10b show the relative frequency (probability distribution) of the polygon area; only the data from Sm-site are plotted separately in an inset, as the typical polygon area at Sm-site is significantly smaller than those at the other sites. All the probability distribution curves show upwardly convex curves with a distinct peak. However, the positions and widths of the peaks vary greatly depending on the locality. For example, the O-site curve shows a peak at ca. 2.5 m², implying that the column diameter is roughly equal to the length with arms extended sideways. In contrast, the curve of Sm-site shows a peak at ca. 2.1 × 10⁻³ m², indicating that the polygon area is approximately the same as the area of an adult's palm. Figure 10c shows the distribution of a normalized polygon area, that is, the measured area A divided by the mean value ${\mu }_{\text{Area}}$. All data points appear to lie on a standard curve, with the exceptions of the Y-site data's peak position, which is marginally shifted to the left, and Sz-site data's peak width, which is narrow. It should be noted that in Figures 10a–10c, the typical polygon area decreases in the order O-site (andesite), Y-site (basalt–basaltic andesite, or alkaline basalt and trachybasalt), Sm-site (rhyolite), and Sz-site (andesite). This suggests that the specific cross-sectional area of columnar joints is not strongly influenced by SiO₂ content, which is known to be rich for rhyolite and poor for basalt.

Figure 10d shows the relative frequency of the number of vertices. This plot demonstrates the predominance of hexagons and pentagons at all the localities. More precisely, the O, Y, and particularly the Sz sites are dominated by hexagons, while the Sm-site predominantly yields pentagons. The probabilities of observing polygons with squares, heptagons, and octagons are significantly small, less than 0.2 for each. The predominance of hexagons and pentagons at the localities is consistent with previous field observations at other study sites (Beard, 1959; Phillips et al., 2013; Walker, 1993). Their predominance has been explained theoretically using the principles of fracture mechanics and the finite element method (Hofmann et al., 2015). Specifically, the dominance of pentagons only at the Sm-site would be attributed to the presence of many T-junctions, as mentioned in Section 2 because the number of vertices of a polygon that involves nearly vertical joint angles tend to be smaller than that of polygons that consist mainly of obtuse joint angles. This point will be revisited in Section 4.2.

Figure 10e shows the relative frequency of joint angles at crack intersections. We determined a peak within the range of 1.8–2.3 radians for all the data. This result is consistent with that given by Figure 10d, confirming that hexagons and pentagons are predominant. In addition, the peak width for the Sz-site curve is as narrow as the peak width of the Sz-curve for the cross-sectional area in Figure 10c, indicating a high degree of regularity in the geometry of crack patterns at the Sz-site.

Coefficients of variation (CV) values, defined as the standard deviation $\sigma$ divided by the mean value $\mu$, are summarized in Figure 11 for all 12 single-peak curves plotted in Figure 10. The horizontal axis shows the columnar joint location in increasing order of the typical cross-sectional area. Note that Figure 11 includes an estimate of CV values for columnar jointing at Giant's Causeway (GC; Northern Ireland), which is one of the most famous columnar joints in the world. It is located along the north coast of Northern Ireland and mainly composed of basalt. The CV of GC was estimated from the line drawing in O'Reilly's paper, O'Reilly (1879), to which we applied the image processing technique described in Section 3.

It is evident from Figure 11 that spatial variation in all three geometric attributes is most suppressed in GC. The plots also show an increasing trend in the CV values with increasing typical cross-sectional area, except for the Sm-site. The large fluctuations in the three geometric properties are due to the difference in the degree of regularity of the crack network. According to the coupled conductive and convective heat transfer model (Hardee, 1980; Budkewitsch & Robin, 1994), the deeper the cracks penetrate inside of the lava, the more the crack network becomes a regularly ordered columnar structure with limited diameter variation as described in Figure 1. Moreover, T-junction gradually changes to Y-junction as the crack advances into the inner lava. As a consequence, the regular crack network is expected to have a more uniform sectional area and shows a hexagonal pattern caused by Y-junction. Comparing the results of the five sites in Figure 11, the fluctuation of CV values in Sz-site and GC is less than that of the other three sites for all geometric attributes. This means that the crack pattern in Sz-site and GC are well ordered with high maturity. Since the Sz-site was formed within a sill (Izu Peninsula Geopark, 2021), it might be cooled more slowly under slight temperature difference between the lava margin (where is cooled at first) and internal lava than the other sites. Due to this slow cooling process, a regular crack network may have been formed at Sz-site.

Another notable fact is that although the lithologic character of Y-site and GC is commonly classified as basalt, the results of CV value in Y-site and GC show different degree of fluctuation. This fact again implies that lithologic character may not affect the degree of variation in columnar joints' geometric attributes, as was observed for the specific cross-sectional area values.

Correlation Between Geometric Attributes

Our systematic studies on polygonal crack networks enable us to examine the interrelationships between the geometric attributes. Figure 12 presents the correlation diagram between the cross-sectional area and the number of vertices for each locality. Data points are widely distributed in the vertical direction from 3 to 11, while the number of extreme polygons with more than eight or less than four vertices is significantly small. An essential deduction from Figure 12 is a positive trend, which signifies that polygons of greater area tend to have a higher number of vertices, independent of localities (i.e., spatial scale or lithologic composition). The reason for the correlation between the cross-sectional area and the number of vertices can be understand intuitively. In a typical polygon pattern observed in a columnar joint, each vertex (i.e., the junction) is shared by three adjacent polygons, and each side (i.e., the crack) is shared by two adjacent polygons. Since the sum of the three internal angles surrounding a vertex is equal to 360°, each of the three internal angles is most likely to have a value close to 120°, which equals to the internal angle at a vertex of a regular hexagon (Figure 13a); this causes hexagons to occur most frequently. However, for a polygon with more than six sides, the average of the internal angles at a vertex exceeds 120°. Therefore, the two other adjacent polygons that share the same vertex as this polygon with more than six sides will have an internal angle of less than 120° (Figure 13b). This fact leads us to the following deduction; the more vertices a polygon contains, the larger the area of that polygon, provided the edge length is nearly uniform within the whole pattern. This deduction is consistent in quality with the experimental data shown in Figure 12.

Deviation from an Ideal Polygonal Pattern

Using our obtained statistical data, we conduct an additional analysis to elucidate the geometric deviation of the naturally occurring fractured pattern from an ideal regular polygonal pattern. Our method of analyzing the regularity of polygonal shapes is based on Vasseur's approach, Vasseur and Wadsworth (2019), in which the relationship between the averaged side-length b and the area A is examined. For example, provided the polygons are regular hexagons, all the data points lie on the single curve of A = 2.60b² in the A–b plane, regardless of the polygon's length scales. The relationship between regular polygons with n sides can be proved by elementary geometry and is described by the following: [Image Omitted. See PDF]

Approximate values of the coefficient C_n for n = 4, 5, 6, and 7 are summarized in Table 2. Therefore, by examining the deviation of data points from the ideal curve given by Equation 1, we can evaluate the geometric deviation of naturally occurring fractured patterns from regular ones.

Table 2 Values of C_n and ${\alpha }_n$ for n = 4, 5, 6, 7 to Three Significant Digits

n	4	5	6	7
Cn	1.00	1.72	2.60	3.63
${\alpha }_n$	0.815	0.853	0.865	0.840

Figure 14a shows our four-site measurement data in the A–b plane with a double logarithm plot. The solid slanted lines correspond to the relation of Equation 1 with various n values. Data points are distributed on or below the solid line for each polygon because a regular polygon has the maximum area when the total side lengths are fixed. To quantify the geometric deviation of natural n-gons from regular ones, we draw a dashed slanted line parallel to the solid line such that 95% of the data points are included between the two lines. Here, the outlier points (5% of the whole for each n) correspond to polygons with largely distorted shapes far from regular polygons, as shown in Figure 14c. Note that each dashed line can be represented by $A^{\prime }={\alpha }_n{C}_n{b}^2$ with an n-dependent constant ${\alpha }_n$. Consequently, the magnitude of ${\alpha }_n$ quantifies the regularity of the naturally occurring n-gons. ${\alpha }_n$ equals to 1 when all n-gons contained in the natural pattern are entirely regular. Conversely, ${\alpha }_n$ may be close to 0 when a subset of n-gons is exceptionally distorted.

The obtained values of ${\alpha }_n$ are listed in Table 2. They are approximately 0.85, regardless of the polygon type. This result implies another common rule that the columnar joints follow. Namely, the geometric distortion of polygons that appear on the columnar joints' outcrop is suppressed to ∼0.85 when ${\alpha }_n$ is defined as above. This common rule holds over a broad spatial scale ranging from a centimeter-sized polygonal pattern (at Sm-site) to a meter-sized pattern (at Y-site).

Discussions

Numerous table-top analog experiments have successfully reproduced the polygonal prismatic shape of columnar joints. Most of these experiments are based on desiccation cracks of starch-water mixtures (Akiba et al., 2017; Ellis & Blenkinsop, 2019; Ma et al., 2019; Müller, 1998a, 1998b) and calcium carbonate-water mixtures (Akiba & Shima, 2019; Nakahara & Matsuo, 2006). The mixtures consist of solid grains and pore spaces filled with either liquid or air bubbles. On drying, the liquid content evaporates, causing volumetric shrinkage of the mixture. Subsequently, cracks occur at the surface and develop toward the interior of the mixture. As a result, these materials exhibit a regular array of tiny polygonal prisms reminiscent of columnar joints. Additionally, stearic acid (Christensen et al., 2016) and metachalk (Weinberger & Burg, 2019) are also known to form prismatic patterns when they are cooled. Despite the differences in chemical composition and scale length, the apparent crack patterns are relatively similar to those of columnar joints. This similarity arises from physics' general laws that govern the volumetric shrinkage of both cooling lava and drying mixtures (Müller, 1998a). The next paragraph provides the rationale for supporting the similarity of volumetric shrinkage between the two systems.

Suppose that a cooling layer of lava is subjected to the constrained conditions of no horizontal displacement. Then the thermal contraction stresses in the horizontal direction read [Image Omitted. See PDF]with $\mu$ being the shear modulus; the vertical component of the stress, p_z, vanishes under the plane stress approximation. The volume contraction term in Equation 2 is given by [Image Omitted. See PDF]where $\sigma$ is Poisson's ratio and [Image Omitted. See PDF]is the thermal volume contraction under unconstrained conditions; $\alpha$ is the volume expansion coefficient, and T(z, t) is the spatio-temporal temperature profile (Turcotte & Schubert, 1982). Equation 4 means that the stress occurs only after cooling below T₀, where T₀ was estimated as ca. 1000°C, being close to the solidus temperature of lava. The temperature profile evolves over time according to the diffusion equation given by [Image Omitted. See PDF]where the thermal diffusivity ${\kappa }_{\text{th}}$ of lava is about 10⁻⁶ m²/s. Using typical material parameters of lava: $\mu$ = 33 GPa, $\sigma$ = 0.25, and $\alpha$ = 10⁻⁵ K⁻¹, and tensile strengths of 10–30 MPa (Lockner, 1995), the temperature differences T(z, t) − T₀ was estimated as −30°C–−100°C for the crack front. Once the tensile stress driven by the temperature difference exceeds the fracture strength of solidified lava, the new crack appears at the existing crack front.

Quite interestingly, the parallel argument as above holds true for dried starch slurry. The horizontal tensile stresses in a drying starch layer are also described by Equation 1, and the volume contraction is [Image Omitted. See PDF]where ${\left(\mathrm{\Delta }V/V\right)}_d$ is the unconstrained volume contraction due to desiccation. This quantity can be expressed by the water concentration C = V_w/V: [Image Omitted. See PDF]where V_w is the water volume in the total volume V, and the subscript 0 denotes the values prior to desiccation. Using the approximation of $V_0/V\approx 1$, we have [Image Omitted. See PDF]and the obtained result of (8) is analogous to Equation 4. Hence the water concentration profile C(z, t) plays the same role for desiccation as the temperature profile T(z, t) for cooling lava. Moreover, the water concentration profile also obeys the diffusion law: [Image Omitted. See PDF]where the hydraulic diffusivity ${\kappa }_h$ of starch is 10⁻⁸ m²/s. In summary, the contraction mechanisms of cooled lava and that of dried slurries are essentially similar.

The analogy between cooling lava and drying mixtures can also be applied to the relationship between the crack spacing (or section area) and the material thickness. In the drying particulate suspensions, the crack spacing, $\lambda$, is proportional to the power of material thickness, h, as $\lambda \propto h^{\beta }$ (Smith & Sharp, 2011). The power law also holds for the polygonal area, A, as $A\propto h^{4/3}$ (Flores, 2017). For many different materials, the power laws were found to hold for any value of h. A notable exception is the case of dried cornstarch slurry, in which the polygonal area converges to a constant when the sample thickness exceeds a certain critical value (Akiba et al., 2017; Ma et al., 2019). This converging behavior of dried cornstarch slurry is similar to that of cooled lava, in which the column width is believed to be constant when the crack propagation exceeds a certain depth from the lava flow margin (DeGraff & Aydin, 1993; Goehring & Morris, 2008). In fact, for columnar joints, a characteristic depth below which the stria size (which is proportional to the column's cross-sectional area) converges to a constant was theoretically predicted as several meters; this value was obtained from the estimation of the thickness of the heat transfer layer (Goehring & Morris, 2008). In view of the material thickness dependence, therefore, both the cooled lava and dried slurry show the similar converging behavior in the evolution of the crack spacing (and also the cross-sectional area).

Results of the starch-based analog experiment (Toramaru & Matsumoto, 2004) provide the following two conclusions. First, the average cross-sectional area of starch columns is inversely proportional to the drying rate. Therefore, the average area decreases under fast drying conditions and increases at slow drying rates. Second, the preferred column shape also depends on the drying rate. Pentagons occur most frequently under fast drying conditions, while hexagons occur most often at slow drying rates. Based on the two conclusions, we can say that the three physical factors (i.e., drying rate, average area, and polygon type) are interdependent in the starch columns (see Figure 15a). Furthermore, it is naturally inferred that a similar interdependence exists in columnar joints, whereby the cooling rate of lava replaces the drying rate of moist materials. Some of this speculation has already been confirmed. It was experimentally observed that fracture advance speed and fracture spacing in columnar joints are inversely related (Goehring et al., 2009). This indicates that large (or small) column diameters are caused by slower (or rapid) cooling rates, consistent with conjecture from crystal texture analysis (Long & Wood, 1986). In addition, the similarity in the contraction mechanism mentioned earlier implies that the correlation between the contraction rate and polygon shape, experimentally observed in dried slurries, will hold true for the cooled lava, too.

Before the completion of the present study, however, it was unclear whether the number of vertices in a cross-sectional polygon is related to the column diameter or cooling rate. Against this background, we scrutinized whether the number of vertices in columnar joints is positively correlated with column diameter size, as demonstrated in Figure 12c. This finding suggests that a three-cornered relationship exists in columnar joints (see Figure 15b), similar to the starch-based analog experiments. More concretely, we established that columnar joints' cross-sectional shape shows the following tendency: pentagonal columns with small diameters tend to occur frequently at rapid cooling rates, while hexagonal columns with large diameters are more common at slower cooling rates. Field measurements or numerical simulations are required for verification.

Conclusions

This study demonstrated the statistics and intercorrelation of geometric attributes that characterize crack patterns observed in Japan's columnar joints. Drone photography and computer-based image analysis of the polygonal crack patterns provided highly accurate data on the area, vertices, and inner angles. Data derived from four different localities showed common features in the occurrence frequency of the latter two attributes and the normalized polygon area. However, the data revealed a large discrepancy in the typical polygon area. Moreover, we established a positive correlation between the cross-sectional area and the number of vertices, irrespective of lithologic character and spatial length scale of the columnar joints. Our findings on columnar jointing agree with previous studies using numerical simulations and starch-based analog experiments, suggesting the existence of universal geometric features regarding polygonal crack patterns.

Acknowledgments

Y. Akiba appreciates the financial support from the Sasakawa Scientific Research Grant from the Japan Science Society. This research was also supported by the Japanese Society for the Promotion of Science Grants-in-Aid for Scientific Research (JSPS KAKENHI, grant nos. 18H03818, 19H05359, 19K03766, and 20J10344). Fruitful discussions with Motohiro Sato are gratefully acknowledged. The authors' great thanks are also extended to Hideo Hoshizumi, Isoji Miyagi, Kei Kurita, and Ichiro Kumagai for their productive discussions and valuable suggestions. Finally, they would like to express our warm appreciation to all the kind staff in the following public agencies for their generous cooperation during the surveys and data processing: Hagi Geopark Promotion Council (Y-site), Shimoda city office (Sz-site), Izumo-shi government office (Sm-site), and Kumejima-cho public office (O-site).

Data Availability Statement

The datasets used in this study will be stored on Figshare, and their copies are temporarily uploaded as supporting information for review purposes (Supporting Information S1.pdf, Table S1.csv).