Understanding the spatio-temporal evolution of fractures in pillow basalt

We investigated the origin and spatio-temporal evolution of cooling fractures in pillow basalt which undergo thermal contraction after their eruption in an aqueous environment. Through a computer-based simulation using Fourier transformation, the thermo elastic stress displacement profiles within individual pillow units are determined. The scaled model (pillow diameter - 1 meter) generated radial, linear fractures perpendicular to pillow margin and irregular discrete flaws in the pillow interior like the ones observed in natural examples. Radial linear fractures of 3–5 centimetre in length have been measured in pillows of average one-metre diameter from the Maradihalli region, in the Chitradurga Schist Belt, India. An estimated time of 94–118 minutes was required to get radial fractures of similar length in the simulation. Our model efficiently replicated the generation and distribution of thermal fractures and allowed an estimation of cooling time for the peripheral glassy zone but has limitations in deciphering the formation of fracture networks in progressively crystalline inner zone of pillows.


INTRODUCTION
'Pillow basalt' is a term commonly associated with lava flows of basaltic or similar compositions made of bulbous, ellipsoidal or tubular lobe masses about the form and size of pillows.Some researchers have considered pillow basalt to be composed of discrete, independent, sack like or ellipsoidal masses (Johnston, 1969;Macdonald, 1972;Snyder and Fraser, 1963) and others regarded pillow basalt as a tangled mass of cylindrical, interconnected flow lobes (e.g.Vaugnat and Pustaszeri, 1965).The pillows are considered to form when fluid lava cools in contact with water, either when it erupts directly into water or when it flows into a body of water.(Snyder and Fraser 1963).Moore (1975) suggested a mechanism of formation of such structures based on real time documentation of the process by still and motion picture photographs of lava from new Mauna Ulu vent of the Kilauea volcano in Hawaii.The study revealed that the lava flow is carried out by a tube of lava, whose outer carapace is quickly chilled and forms a thin crust when it comes in contact with water but inside it has molten lava.The diameter of individual pillows ranges from several centimeters to meters.Most of the pillows commonly exhibit radial fractures inside them which result from thermal contraction with progressive cooling (Wells et al., 1979).Considerable amount of research has been conducted on the mineralogy and geochemistry of pillow basalt and their We investigated the origin and spatio-temporal evolution of cooling fractures in pillow basalt which undergo thermal contraction after their eruption in an aqueous environment.Through a computer-based simulation using Fourier transformation, the thermo elastic stress displacement profiles within individual pillow units are determined.The scaled model (pillow diameter -1 meter) generated radial, linear fractures perpendicular to pillow margin and irregular discrete flaws in the pillow interior like the ones observed in natural examples.Radial linear fractures of 3-5 centimetre in length have been measured in pillows of average one-metre diameter from the Maradihalli region, in the Chitradurga Schist Belt, India.An estimated time of 94-118 minutes was required to get radial fractures of similar length in the simulation.Our model efficiently replicated the generation and distribution of thermal fractures and allowed an estimation of cooling time for the peripheral glassy zone but has limitations in deciphering the formation of fracture networks in progressively crystalline inner zone of pillows.environment of formation (Bishop et al., 2002;Cruz et al., 1989;Ramsay and Moore, 1985;Yazar, 2018;Zhang and Zhu, 2018).Borradaile and Poulsen (1981) and Borradaile (1982) documented tectonic deformation of pillows and their implication in understanding the true bedding and way-up direction.However, few studies have addressed the internal structure and fractures related to thermal stresses within individual pillow units which have an important bearing on their genesis and thermal evolution.Bergbauer et al. (1998) have computed the thermo-elastic stresses from the temperature profile (conductive heat transfer) of 2D domains using Fourier transformation.The relation of the elastic stresses to the thermal potential is described by Timoshenko and Goodier (1951).However, the calculations are not generalized for anisotropic systems.In light of this, the paper attempts to trace the evolution of thermal fractures within individual pillow units from Maradihalli, Chitradurga, Karnataka, India through a 2D thermo-mechanical modeling.The generalized mathematical framework can be adapted for any irregular shaped structural objects and anisotropic material properties.
Craton is known to preserve ancient Archean continental crust >3.0Ga (Jayananda et al., 2006).The Eastern and Western blocks of the Dharwar craton (EDC and WDC) were juxtaposed against one another between 2.75-2.51Ga,along the Chitradurga Shear Zone (CSZ; see Figure 1, e.g.Chadwick et al., 2003;Jayananda et al., 2006;Naqvi and Rogers, 1987).The volcano-sedimentary sequences of the CSB in the eastern part of WDC are metamorphosed up to green schist to lower amphibolite facies (Chadwick et al., 1989;Jayananda et al., 2006;Mondal, 2018;Mondal and Mamtani, 2014).The pillow basalt from the Chitradurga Group (Chakrabarti et al., 2006;Roy et al., 2020) was formed ~2.8Ga ago when hot basaltic magma extruded from submarine vents and cooled rapidly on the ocean floor.Continuous eruption led to the formation of bulbous structures (Fig. 2B, C), generating interconnected pillow basalt units, whose occurrence within a back-arc basin (Bhowmick and Mondal, 2021;Manikyamba et al., 2021) is indicative of uninterrupted lava flow.A cross-section of an individual pillow shows distinct colour contrast separating the chilled margin, the glassy zone and the relatively crystalline inner zone.Progressive crystallization at slower cooling rates toward the interior reduces the amount of glass and makes the centre of the pillow more crystalline.Vesicles are mostly concentrated along the glassy part of individual pillows.A plethora of radial fractures of varying lengths exist along the boundary and within the body of the pillows (Figs.2C; I Appendix).These fractures mostly terminate at high angles to the pillow boundary which is indicative of their thermal origin.The fractures are closely spaced and mostly within the glassy zone.Within this zone, the fractures are linear, and from a few millimetres to a few centimetres in length (Fig. 2C inset).Towards the interior of the pillow, the fractures become irregular.There, fractures are relatively widely spaced forming a network of polygons.The pillow basalt also shows prominent vesicle pipes perpendicular to the glassy selvage.When a glassy selvage has developed on a pillow, the slower cooling within the pillow causes crystallisation, which also entails the release of any gas the magma may have contained.Gas bubbles begin to form as a result, and when the crystallisation front proceeds inward, the exsolved gas is joined to the initially formed bubbles to create vesicle pipes (Philpotts and Lewis, 1987).

METHODOLOGY
The heat transfer from the centre of the pillow towards the periphery is by conduction due to the thermal gradient.Initially as the pillow has emerged due to volcanic activity, it is at a high temperature.With time, the solid pillow cools down as the heat energy slowly diffuses (lost by transfer) to the surrounding water.Since the conductive resistance inside the pillow is higher than on the solid surface-fluid interface, the spatial temperature distribution inside the solid pillow is important.Here we study the temperature dynamics inside the solid pillow and the surrounding water body.The shape of the pillow is considered an ellipsoid, having the base diameter of 1m and height of 0.6m.The schematic of the problem geometry is shown in Figure 3A.
The spatio-temporal temperature (T) evolution in the pillow is described by the conductive heat transfer equation (Bird et al., 2015) including the latent heat of crystallization (Q), (1) where ρ is density, C p is the specific heat capacity, k is the thermal conductivity.The liquidus and the solidus temperature of the molten (basaltic) magma is 1200ºC and 1030ºC (Liu and Lowell., 2011).Depending on the molten phase fraction (θ), the physical properties of the mixture phase can be correlated (Lienhard, 2005) where 1 and 2 represent the molten phase and the solid phase.Further it is assumed that since the melt viscosity is so large, the flow dynamics of the molten phase is insignificant within the small spatial dimension and the lesser cooling time scale.Thus, the buoyancy driven flow of the molten magma during the cooling phase in the pillow is not considered.The heat flux on the surface of the pillow, in contact with the surrounding (cold) water can be described by the nucleate boiling flux (Rohsenow correlation; Som, 2008) as, (6) where n is the outward unit surface normal, the subscripts L and V represent the liquid and vapour phase of water, Q LV is the enthalpy (latent heat) of vaporization of liquid water, σ is the liquid-vapour interfacial tension, T sat is the water saturation (boiling) temperature, Pr L is the Prandtl number, which varies from 2-10 depending on the water temperature.The parameter C sf is the surface-fluid factor, and we have considered a typical value of 0.006 in this case (Som, 2008).The base of the pillow is considered at far field condition, ∇T=0.
The heat transfer in the surrounding water is accompanied by the buoyancy driven flow, which can  Once the temperature distribution is obtained, the stress displacements can be computed using equation ( 14).
The conductive time scale (τ) of the system is L c 2 /α, where L c is the characteristic length scale (base radius of the pillow in this case).This suggests that the temperature dynamics is relevant for t~10τ, beyond which it attains a steady state (cools down).However, accounting for the latent of crystallization, the cooling duration is increased by 20-30 times, which suggests the cooling period of the pillows (Nichols et al., 2009), and the active time window of analysing the thermal stress formation.
Modelling the water domain, the dimension of the surrounding box is such that it is large enough to allow negligible heat transfer across the domain boundaries (no thermal gradient, far field condition).We have solved the coupled system of equations numerically using the finite element method with the help of COMSOL Multiphysics package v6.0.
The discretization of the variation form of the governing equations follows the standard Galerkin formalism.To save computational cost, the Jacobian matrix is updated once in several iterations.Within each Newton iteration, the sparse linear system is solved by preconditioned Krylov methods such as the generalized minimum residual (GMRES) method and the bi-conjugate gradient stabilized (BCGSTAB) method.The weak boundary condition is solved using an inexact Newton method with backtracking for enhanced convergence and stability.
The fully coupled solver is used with a relative tolerance of 0.001.The modelling domain is meshed with free triangular elements (first shape order) and around 600,000 degrees of freedom are solved.Boundary layer meshing of five layers is constructed near the solid-liquid interface (pillow surface).Mesh independence was tested at various mesh sizes from 0.0008 to 0.1m.It was found that the total heat flux from the surface of the pillow is unchanged (variations less than 10 -4 ), for mesh size less than 0.0015.Hence, the mesh size (maximum) of 0.015 was considered adequate for the current numerical simulations.Evolution of the stress displacement profiles are observed in xz cut planes (Fig. 3B) at different time instants.

Stability of the governing equations
There are two major set of equations which are numerically solved in this work.The first one is the thermal diffusion (conduction) inside the pillow as described by equation 1 and the second is the natural convection driven flow surrounding the pillow, described by equations 7-9.
Here we present the normal mode stability analysis of the physical system of equations.The general idea of linear stability analysis is to introduce an infinitesimal perturbation to a background state and examine whether its amplitude grows or decays with time.Since the perturbations are small, the equations can be linearized which aids the mathematical analysis.In the method of normal modes, the perturbations are assumed to be sinusoidal.The differential equations 18 and 19 are linear and the coefficients are spatially and temporally invariant.Thus, the normal mode ansatz for and can be described (Nayfeh, 2008) as, where a is the wave number and σ is the growth rate of the perturbation.If σ>0, the disturbances will grow in time and the system will be unstable.Similarly, for σ<0, the perturbations will subside in time leading to stable solution.For σ=0, the marginally stable case, is the limiting state distinguishing between the stability and instability regimes.Thus, substituting equations 22 and 23 into equations 20 and 21 leads to This is an eigen value problem with being the eigen function.The eigen values for a given a will give Ra for a (limiting) stable solution of the equations of motion.The solution of the eigen value problem is dependent on the type of boundary conditions.In the case of free-free boundary, the value of Ra(a) for marginal stability of the n th node of wavenumber a (Drazin, 2002) is, (27) This gives the demarcation between the growing and decaying modes in the plane of Ra and a, such that the minimum value of the Ra on the curve is its critical value.The minimum can be obtained for n=1 (least stable vertical mode) and then the minimum for all real horizontal wavenumbers a. (28) Similarly for the free-rigid and rigid-rigid boundary conditions, the critical Ra is 1101 and 1708 respectively, for unconditional stability of all possible wavenumbers.The Ra used in this work is much less than the critical value, thus there is no possibility of convective instability and hence, the unphysical evolution of the flow and temperature fields in the numerical simulation are not expected.The courant number (Pletcher et al. 2012) for the computation is maintained at less than 1, ensuring numerical stability of the calculations with convective transport.

Conductive heat transfer inside the pillow -stability analysis
In the case of uniform thermal conductivity, equation 1 can be recasted as t a , 2 1 . 8 , 1 -1 0 , I ( 2 0 2 3  Spatio-temporal evolution of fractures in pillow basalt 7 where is the thermal diffusivity and is the source term accounting for the phase transition.We are showing here the analysis for a simplified 1D scenario.However, the results are also applicable in the case of 2D and 3D settings.Adding a small perturbation to the base State (T 0 ), the perturbation equation to the first order is ( for). (30) The perturbation variable can be once again expressed as the normal mode ansatz (refer to Eq. 23) such that (31) Substituting this into equation (30) leads to (32) Since the source term is related to the solidification cooling of the molten pillow (crystallization), f~-T.Thus, f'<0 and consequently σ<0.Thus, for all possible wavenumbers m, the system is inherently stable and all disturbances are suppressed with time.Numerically such systems do not exhibit any unphysical evolution of temperature fields and are not sensitive to numerical settings.The von-Neumann stability condition is also ensured in the computational settings for diffusive systems in this case.

RESULTS
The thermo-physical properties of the pillow basalt are considered based on the data reported in the work by Liu and Lowell (2011).The thermo-physical properties of water are taken from Perry's handbook (Perry, 1950).
Since the computational results are in the 2D domain, we have chosen one cut plane view to present the results.Figure 4 depicts thermal stress displacements (flaws) as generated in the simulated models.The regular linear radial flaws are observed, terminating at high angle to the pillow margin which show systematic increase in length with progressive cooling.
As the hot lava inside the pillow quenches, it contracts, generating a number of small, regular fractures radially along the chilled margin (Fig. 4).With progressive cooling the radial fractures along the margin propagate linearly towards the centre of the pillow and their length increases up to a certain limit.At various time intervals, the xz cut plane is used to illustrate the evolution of the stress displacement profiles.The xz planes are displayed in time instants t= 0 minutes to 117.99 minutes in 4b(a) to 4b(i), respectively.

DISCUSSION
The outer surface of a natural pillow forms a very thin crust as soon as it comes in contact with the surrounding water whose ambient temperature is drastically lower than lava temperature.This high temperature gradient leads to very rapid cooling, thus inhibiting crystal growth, resulting in the formation of a glassy outer rim, immediately adjacent to the chilled margin.As soon as the glassy rim develops, it causes an insulating effect around the inner portion of the pillow which is still filled with hot lava.Thus, the temperature gradient in the interior of the pillow falls, leading to slower cooling rate, resulting in the formation of crystals.It must be emphasized here, that there are two distinct timescales in the problem, one is related to the cooling of the pillow and the other is the thermal fracture (particularly near the peripheral region).From our calculations, it is evident that the fracture formation due to the thermal expansion (2-10hours) is much faster compared to the cooling due to the thermal diffusion (5-10days).Since the volume shrinkage is primarily associated to the cooling, it is unlikely that there is significant volume shrinkage occurring during the development of the thermal stress induced fractures.
Field photographs show a distinct colour contrast separating the outer glassy zone from the inner more crystalline zone (Fig. 2B).It may be noted that the regular and linear radial fractures are mostly limited to the outer glassy zone.Some of the radial fractures extend across the glassy zone, becomes irregular in orientation and tend to get inter-connected forming a network of polygons in the internal portion of individual pillows (Fig. 2C).The irregularity of the fractures beyond the glassy zone can be attributed to a control of increasing crystal fraction towards the interior of the pillows, which inhibits linear propagation of fractures.
The result of the simulated model closely approximates the natural examples with one major exception.The linear radial fractures along the pillow margin and the irregular discrete flaws within the pillow interior generated in the simulation always remain in two distinct domains and they never get interconnected.This contrast of fracture pattern can be attributed to the fact that the simulated model only considers the difference of heat diffusivity/conductivity between the inner and outer parts of the pillow but does not incorporate the control of textural difference between the outer glassy and inner crystalline parts of the pillow.pillows of average one-meter diameter which are generated at 90 to 100 minutes time instants of the simulated model (Fig. 5).Therefore, an estimation of cooling time for the outer glassy zone of natural pillows can be inferred provided the assumed thermo-physical constants remain uniform over the specified period of time.It is worth mentioning that the pillows observed in the Maradihalli region are not in their original shape, since they underwent multiple deformations.Therefore, the orientation of fractures within them as observed now will not properly reflect the simulated fractures in an ideal pillow body of regular undeformed shape.

CONCLUSIONS
The present study aims at understanding the formation of contractional fractures and their spatio-temporal evolution in pillow basalt units through field investigation and numerical simulation.The major conclusions of the study are as follows: i) The simulation effectively reproduces the linear radial fractures along the pillow margin and the discrete irregular fractures in the pillow interior as observed in natural examples.However, the model is limited in providing a suitable mechanism for the formation of fracture networks in the progressively crystalline inner zone of the pillows since it does not incorporate the presence of vesicles and textural inhomogeneity inside the pillows that develop with time during cooling.
ii) Time required for the generation of contractional fractures within individual pillows can be approximated by comparing the lengths of marginal linear fractures in natural pillows with the ones simulated in the model.
Fig. c FIGURE 3. A) Domain model of a pillow and surrounding water enclosure.B) View of the plane slice (xz) cut at the pillow center.
be described by the coupled momentum and energy conservation equations.For low Reynold number flows, the inertia is negligible, and thus the conservative form of the momentum equation(Fox et al., 2020) is, where v is the fluid flow velocity field.Ignoring fluid heat dissipation, the convective-diffusive heat transport equation can be described(Incropera and DeWitt, 2002) as, (9)We use the Boussinesq approximation to account for density variations(Incropera and DeWitt, 2002),(10)    where β L is the thermal expansion coefficient of water.We consider no slip (flow) boundary condition on the pillow surface, and the seabed.The boundary surface opposite to the seabed is considered an open boundary.All the remaining boundary surfaces are considered periodic condition.Initially, the pillow basalt and the water temperature are considered 1200ºC and 4ºC (cold sea water).The stress displacement u is dependent on the plane stress displacement potential (ψ)(Bower, 2009) as,   (11)    Considering a linearly elastic material (thermal strain ∈~-βT), the thermo-elastic stress displacements due to the thermal gradients in the solid fraction of the pillow, can be obtained based on the theoretical derivations byTimshenko  and Goodier (1951),   (12)    where β is the thermal expansion co-efficient of the solid and v is the Poisson ratio of the material.Taking the temporal derivative on both sides to equation (12), leads to (13)In the case of single-phase material (after crystallization) for uniform thermal diffusivity , we get,(14) D O I : 1 0 . 1 3 4 4 / G e o l o g i c a A c t a 2 0 2 3 .2 1 .8 Spatio-temporal evolution of fractures in pillow basalt 6 Natural convection flow of water outside the pillowlinear stability analysis On non-dimensionalization the fluid flow and the energy equationderivative, Pr=v L /k L is the Prandtl number, is the Rayleigh number.Here L c is the characteristic length scale which can be considered to be the pillow diameter.The other variables are nondimensionalised with the appropriate scales.Considering the Boussinesq approximation (Eq.10), the continuity equation can be simplified to (17) Taking the curl of equation 15 and considering the vorticity ω=∇×v, we get (Drazin and Reid 2004) number of unknown variables, we express the flow velocity components (in 2D) as the derivatives of the scalar stream function and .The vorticity takes the form ω=-∇ 2 ψ and thus equations 18 and 19 can be recasted as case of σ→0, and eliminating from the above set of equations yield,(26) ) D O I : 1 0 . 1 3 4 4 / G e o l o g i c a A c t a 2 0 2 3 .2 1 .8 S .M o n d a l e t a l .

FIGURE 4 .
FIGURE 4. Evolution of the stress displacement profiles shown in the xz cut-planes at progressive time instants.The results shown are for the thermal stress displacement of a single isolated pillow.

FIGURE 5 .
FIGURE 5.The error bars represent the standard deviation of the averaged displacements originated from the pillow surface (periphery)