Quantitative Evaluation of the Fracturing State of Crystalline Rocks Using Infrared Thermography

The fracturing state of rocks is a fundamental control on their hydro-mechanical properties. It can be quantified in the laboratory by non-destructive geophysical techniques that are hardly applicable in situ, where biased mapping and statistical sampling strategies are usually exploited. We explore the suitabilty of infrared thermography (IRT) to develop a quantitative, physics-based approach to predict rock fracturing starting from laboratory scales and conditions. To this aim, we performed an experimental study on the cooling behaviour of pre-fractured gneiss and mica schist samples, whose 3D fracture networks were reconstructed using Micro-CT and quantified by unbiased fracture abundance measures. We carried out cooling experiments in both controlled (laboratory) and natural (outdoor) environmental conditions and monitored temperature with a thermal camera. We extracted multi-temporal thermograms to reconstruct the spatial patterns and time histories of temperature during cooling. Their synthetic description show statistically significant correlations with fracture abundance measures. More intensely fractured rocks cool at faster rates and outdoor experiments show that differences in thermal response can be detected even in natural environmental conditions. 3D FEM models reproducing laboratory experiments outline the fundamental control of fracture pattern and convective boundary conditions on cooling dynamics. Based on a lumped capacitance approach, we provided a synthetic description of cooling curves in terms of a Curve Shape Parameter, independent on absolute thermal boundary conditions and lithology. This provides a starting point toward the development of a quantitative methodology for the contactless in situ assessment of rock mass fracturing. We use Infrared Thermography (IRT) to quantify the cooling behaviour of rocks with different fracturing states in the laboratory. Rock samples characterized by higher fracture intensity / porosity cool faster than less fractured rocks. We model experimental cooling curves with a lumped capacitance equation and derive a Curve Shape Parameter independent on boundary conditions and lithology. The Curve Shape Parameters shows clear correlations with fracture abundance measures. Results lay the foundations to develop an upscaled methodology to predict the fracturing state of rock masses. We use Infrared Thermography (IRT) to quantify the cooling behaviour of rocks with different fracturing states in the laboratory. Rock samples characterized by higher fracture intensity / porosity cool faster than less fractured rocks. We model experimental cooling curves with a lumped capacitance equation and derive a Curve Shape Parameter independent on boundary conditions and lithology. The Curve Shape Parameters shows clear correlations with fracture abundance measures. Results lay the foundations to develop an upscaled methodology to predict the fracturing state of rock masses.


Introduction
The fracturing state of rocks is a fundamental control on their geomechanical properties including strength, stiffness and permeability (Oda 1985;Hoek and Brown 1997;Hoek and Diederichs 2006;Paterson and Wong 2005;Rutqvist 2015), which in turn strongly affect the hydro-mechanical 1 3 behaviour of fractured media at any scale of consideration Brown 1988, 1997;Rutqvist and Stephansson 2003). Fracturing is also the main indicator of the damage state of rock masses that tracks brittle deformation processes across faults (Brideau et al. 2009;Peacock et al. 2017), around underground excavations (Cai et al. 2004a, b;Martin and Chandler 2004) and in natural and artificial slopes (Eberhardt et al. 2004;Agliardi et al. 2013;Riva et al. 2018). For these reasons, a reliable assessment of the fracturing state of rocks is key to several geological, engineering and geohazard applications. These include the understanding of natural and induced seismicity (Gu et al. 1984;Rawling et al. 2001), the characterization of fractured reservoirs for water, hydrocarbon and geothermal exploitation (Deshowitz and Miller 1995;Gischig et al. 2020) and the assessment of rock mass properties for stability modelling of slopes and underground excavations in rocks (Hoek and Brown 1997;Pine and Harrison 2003;Wyllie and Mah 2004;Brideau et al. 2009).
Current approaches to the quantification of rock fracturing (Dershowitz and Herda 1992) depend on the scale of consideration, ranging from lab samples to fractured media on the in situ scale of geological and engineering problems (Hoek and Brown 1988;Hudson and Harrison 2000;Diederichs 2007). Rock fracturing descriptors vary according to the sampling dimension used (count-1D, linear areal-2D or volumetric-3D) and the topology of the fractures (points, lines, areas, volumes; Dershowitz and Herda 1992).
A complete representation of the fracturing state in a 3D framework, expressed in terms of volumetric fracture intensity (P 32 ) and fracture porosity (P 33 ), can be obtained using non-destructive imaging techniques such as X-Ray Computed Tomography (Keller 1998;Montemagno and Pyrak-Nolte 1999;Ketcham 2006;Agliardi et al. 2014Agliardi et al. , 2017. This high-resolution technique is very powerful for laboratory-scale studies, geophysical techniques are usually applied in situ (e.g. seismic methods, Ground Penetrating Radar; Grasmuek 1996; Bakulin et al. 2000;Dorn et al. 2012). However, geophysical imaging is often unable to support a reliable reconstruction of fracture networks due to lack of sufficient spatial resolution and related difficulties in identifying small fractures that form most of fracture networks and act as key players of fracture connectivity, rock mass permeability and strength.
For this reason, areal or linear (scanline) field sampling strategies (ISRM 1978) are mainly used in situ, yielding fracture abundance descriptors affected by size and orientation biases (e.g. frequency, P 10 , or surface fracture intensity, P 21 ; Zhang and Einstein 1998;Mauldon et al. 2001) or empirical indices of "geomechanical quality" such as the Geological Strength Index (GSI; Hoek and Brown 1997;Marinos and Hoek 2000;Cai et al. 2004a, b;Agliardi et al. 2016). These approaches have different limitations due to the dimensionality of measurements, scale and the effects of heterogeneity, anisotropy, and localization of the fracture patterns that are usually not considered.
InfraRed Thermography (IRT) is an imaging technique that allows the contactless measurement of the temperature of solid surfaces in space and time, based on their emitted radiation in the Infrared bands of the electromagnetic spectrum (Maldague 2001;Vollmer and Möllmann 2018). This technique has been widely applied in engineering to the non-destructive testing of defects, fatigue and damage in industrial materials and engineering structures (Ibarra-Castanedo et al. 2017;Meola et al. 2002;Clark et al. 2003;Aggelis et al. 2010). In the past decades, IRT has been widely applied in geoscience to monitor volcanic and geothermal processes (Neale et al. 2009;Spampinato et al. 2011) and to support geomorphological, hydrological (Price 1980;Moran 2004) and landslide studies (Baroň et al. 2014;Teza et al. 2015;Mineo et al. 2015;Frodella et al. 2017;Guerin et al. 2019). In rock mechanics, IRT has been applied to the study of the thermal response of rocks subjected to stress or water content perturbations (Cai et al. 2020;Zhang et al. 2021), the characterization of porosity , 2022, and the empirical assessment of rock mass quality (Teza et al. 2012;Pappalardo et al. 2016;Mineo and Pappalardo 2016;Grechi et al. 2021). Many applications are qualitative, based on the use of IRT images to identify and map fractures and damage at surface, and only few authors proposed quantitative empirical approaches to relate the thermal response of rock masses to their geomechanical quality related to fracturing . Although these authors showed that the thermal behaviour of rock masses during heating and cooling is related to their fracturing state, a physically based approach to this problem is still lacking.
In this perspective, we designed and carried out laboratory experiments to investigate the quantitative relationships between the fracturing state of rocks and their response to prescribed thermal perturbations. We used Infrared Thermography to monitor the cooling behaviour of pre-fractured cylindrical samples of crystalline rocks (gneiss and mica schist) with a known fracturing state, in both laboratory and outdoor environments. We interpreted our results within a sound theoretical framework to establish quantitative relationships between fracture intensity (and heterogeneity) and the thermal cooling behaviour of fractured rocks, with the perspective to develop a physically based method for the contactless characterization of rock fracturing.

Theoretical Framework
Heat transfer between solid bodies and their surroundings occur by conduction (diffusion), convection and radiation processes (Incropera et al. 1996;Vollmer 2009). Heat diffusion is the result of direct kinetic energy exchanges among particles or quasiparticles within a solid or fluid at rest, depending on the thermal characteristics and size of the object. Such spontaneous heat transfer always occurs following temperature gradients, as described by the second law of thermodynamics. Convection refers to heat transfer at the interface between a solid and a fluid (e.g. air) in motion (Incropera et al. 1996;Vollmer 2009;Bergman et al. 2011), depending on the fluid nature and velocity. Finally, heat transfer by radiation depends on the temperature of a body, its emissivity and the surroundings conditions and atmospheric characteristics, according to the Stefan-Boltzmann law (Vollmer 2009;Vollmer and Möllmann 2018).
While the temporal trend of temperature during heating strongly depends on the thermal forcing (i.e. amount and temporal distribution of heat transferred to the body), the process of cooling (i.e. temperature decrease occurring when no additional heating is applied) follows a relatively simple exponential law, known as the "Newton's law of cooling" (Vollmer 2009). A solid body heated by an artificial (e.g. furnace) or natural (e.g. solar) source tends to balance its thermal state with the surrounding environment by cooling at a rate depending on its thermal properties, chemical composition, geometry, and boundary conditions. A hot body exposed to a cool airstream will transfer heat by conduction across its interior to the surface, and by convection from its surface to the surrounding fluid. At each point in the body, the temperature will decrease until a steady state is achieved (Bergman et al. 2011). In this context, transient heat transfer by diffusion is characterized by the dimensionless Fourier number (Fo) as follows: where: is thermal diffusivity (m 2 s −1 ) t is time (s) L c is a body characteristic length (m) While the role of convective heat transfer is controlled by the dimensionless Biot number: where: h is convection coefficient (Wm −2 K −1 ) k is thermal conductivity (Wm −1 K −1 ) When the Biot number is < 0.1, the internal heat flow is much larger than the heat loss from the surface, and temperature gradients within the body can be considered negligible. In this case, the cooling dynamics of a body with initial temperature T i that cools down to the ambient temperature T ∞ can be described by the simplified "lumped capacitance" equation (Vollmer 2009;Bergman et al. 2011) as follows: where: T i is the initial temperature of the body at t = 0 ( K) T s is the temperature of the body at t ( K) T ∞ is the temperature of the body surroundings ( K) t is time (s) According to the above equations, the cooling rate of a solid body depends on intrinsic material properties, the nature and flow velocity of the surrounding fluid, and a characteristic length (i.e. size or volume to surface area ratio), with smaller bodies cooling faster.
To apply these concepts to a fractured rock system, we hypothesize that the fracturing state affects two main components of the cooling equation as follows: namely he characteristic length (L c ), since higher fracture intensity corresponds to smaller elementary rock blocks; the convection coefficient at interfaces between fractured block, especially when fractures are open. Then, for a material with a given composition and size it should be possible to find a relationship between the fracturing state and the cooling rate.
In this perspective we carried out a laboratory-scale experimental investigation, aimed at verifying that: (a) the cooling rate of a volume of fractured rock depends on its fracturing state, intact block size and fracture heterogeneity; (b) differences between the thermal responses to prescribed thermal perturbations of rocks with different fracturing states can be quantified using Infrared Thermography (IRT); (c) significant quantitative relationships between the cooling behaviour of fractured rocks and fracturing state can be established.

Materials
We studied samples made of two rock types, as follows: gneiss (G, Fig. 1a-g) and mica schist (M, Fig. 1h-n). They are representative of contrasting mineral composition (quartz vs. phyllosilicate-rich), strength and stiffness, and samples have a bulk rock composition that can be considered representative of upper crustal rocks (Rudnick et al. 2003;Agliardi et al. 2017). Both rock types have low intact-rock porosity (< 1.5%), mostly related to disconnected microcracks at the boundaries and across quartz and phyllosilicate grains (Paterson and Wong 2005). The composition, fabric and mechanical properties of the studied rock have been characterized in detailed by Agliardi et al. (2014Agliardi et al. ( , 2017. The studied rock samples were pre-fractured in uniaxial compression by Agliardi et al. (2014Agliardi et al. ( , 2017, who provided a detailed characterization of their fracture patterns, and show a wide range of fracturing degree, making them ideal for our study.

Gneiss Samples (G)
We studied seven cylindrical samples 78 mm in diameter ( Fig. 1). The rock, belonging to the Monte Canale tectonometamorphic unit (Austroalpine Bernina nappe, Central Alps, Italy; Trommsdorf et al. 2005), is an unaltered granodioritic gneiss, made of quartz (30-42 weight percent, wt%), feldspars (35-45 wt%), chlorite (5-13 wt%) and white mica (5-20 wt%) with a total phyllosilicate content less than 20% in weight (Agliardi et al. 2014(Agliardi et al. , 2017. The rock has a gneissic to protomylonitic fabric with a well-developed compositional layering made by alternating microlithons of quartz, plagioclase and K-feldspar (Q-domains), 2-5 mm thick, and mica-rich cleavage domains (M-domains) up to 1.5 mm thick. Quartz and plagioclase form more than 60-80% of the whole rock volume and occur as sub-equigranular aggregates within Q-domains. The main foliation is deformed by tight to gentle harmonic folds inherited from the tectono-metamorphic evolution of the rock, with sub-rounded or rounded hinges and variable wavelength (10-40 mm) and amplitude up to 35 mm (Agliardi et al. 2014). In unconfined compression, the studied gneiss exhibits a medium-low strength and an average Modulus Ratio according to the classification by Deere and Miller (1966). The rock is weaker than other crustal rocks with similar composition, due to the geometric and mechanical controls of anisotropies associated to the folded fabric and exhibits a brittle to semi-brittle behaviour in shallow crustal stress conditions, with complex failure modes including shear fracturing along foliation or along fold axial planes or the development of cm-scale localized brittle shear zones (Agliardi et al. 2014).

Mica Schist Samples (M)
The sample set includes 7 cylindrical samples with a diameter of 83 mm (Fig. 1). The rock belongs to the Edolo Schist Formation, forming most of the metamorphic basement of the central Southern Alps (Spalla et al. 1999), and is an unaltered schist made of white mica (38-45 weight percent, wt%), quartz (25-40 wt%), plagioclase (10-20 wt%) and chlorite (10-15 wt%), with minor amounts of garnet and rutile. The high phyllosilicate content (exceeding 50 wt%) strongly controls rock fabric, that is characterized by a dominant foliation associated with the greenschist facies assemblages. The foliation is mostly continuous at mm-scale, marked by the shape preferred orientation of chlorite and white mica grains, and was formed by decrenulation of preexisting foliations (Spalla et al. 1999;Agliardi et al. 2014). The foliation is deformed by tight to gentle crenulation folds with sub-angular to rounded hinges with wavelength between 10 and 70 mm and amplitude up to 50 mm (Agliardi et al. 2014). Unconfined compression tests revealed lower and less scattered strength characteristics with respect to the gneiss samples, with low to very low strength and high Modulus Ratio Deere and Miller (1966), implying a markedly brittle mechanical behaviour (Agliardi et al. 2014). The rock exhibits fabric-controlled failure modes like those observed for gneiss samples but associated to the development of more spaced and long fractures.

Fracture Network Quantification
The studied samples were pre-fractured by Agliardi et al. (2014) in laboratory unconfined compression tests performed at constant displacement rate using a stiff, servocontrolled loading frame. Brittle failure of both rock types was accompanied by the macroscopic development of shear fractures along the main foliation (planar segments or fold limbs), as well as fractures that cut across the foliation and those that developed parallel to fold axial planes. Fractures formed along foliation or axial planes are generally planar, whereas fractures across foliation have a rough curved or stepped character. More complex failure patterns were also observed, involving the partial or complete development of brittle shear zones exploiting both the foliation and surfaces parallel to fold axial planes (Agliardi et al. 2014).
We used a high-resolution analytical approach, based on X-Ray Computed Tomography, to obtain an unbiased quantification of the fracturing state of the selected rock samples. High-resolution X-Ray CT allows a non-destructive, detailed three-dimensional reconstruction of fracture networks in rock volumes at the laboratory sample scale (Keller 1998;Montemagno and Pyrak-Nolte 1999). The pre-fractured samples were scanned at a resolution of 625 μm using a GE D-600 medical CT scanner (Agliardi et al. 2014). Starting from these image stacks we performed a 3D reconstruction of each sample volume using the software AVIZO™. In these reconstructions, each voxel is attributed a grayscale value proportional to the local X-Ray attenuation, with lowest values corresponding to voids (e.g. open fractures). We identified fracture domains by CT image segmentation using calibrated binary thresholds.
Then, for each sample we selected the fracture domains crossing sample external boundaries (i.e. directly exchanging heat with the outside). Finally, for each sample we quantified fracture surfaces and volumes ( Fig. 2) and, knowing each sample total volume, calculated the corresponding volumetric fracture intensity (P 32 ) and the fracture porosity (P 33; Table 1) according to the definition of Dershowitz and Herda (1992). These descriptors are quantified by sampling fracture domains in a volume space; thus, they provide an unbiased description of fracture abundance (Mauldon et al. 2001;Rogers et al. 2017).

Fig. 2
Example of 3D fracture reconstruction for the sample G16B ( Fig. 1c), performed on the X-Ray CT image stack (voxel size: 625 μm) using Avizo™ software. Fractures are portrayed as rendered solids. 3D reconstruction allowed an unbiased quantification of fracture abundance descriptors P 32 (volumetric fracture intensity) and P 33 (fracture porosity). Their values are summarized in Table 1 Table

Infrared Thermography Experiments
We performed a series of experiments to monitor the thermal response of rock samples, characterized by different fracturing states, to forced heating and subsequent cooling, both in a controlled laboratory environment and outdoor, using Infrared Thermography. Laboratory experiments (Fig. 3a-c) were performed on all the available gneiss (G) and mica schist (M) samples, with the aim of finding relationships between their fracturing states and cooling dynamics under specifically designed controlled conditions. An outdoor experiment ( Fig. 3d) was also carried out on G samples to verify that these relationships can be detected also in uncontrolled environmental conditions, where heating is related to natural solar forcing.

IRT Fundamentals
According to the Planck radiation law, a blackbody (i.e. a body absorbing all incident radiation) emits electromagnetic radiations at different wavelengths depending on its temperature. For earth surface conditions, the peak of radiation falls in the infrared (IR) band of the electromagnetic spectrum. Thus, it is possible to measure the surface temperature of objects based on their emitted radiation in specific IR wavelength ranges, using a technique called Infrared Thermography (Shannon et al. 2005). The thermal radiant exitance, i.e. the energy emitted from a body in the infrared band, is a function of the surface temperature according to the Stefan-Boltzmann law as follows: where: J is thermal radiant exitance (Wm −2 ) T is temperature of the body (K) σ is Stefan-Boltzman constant (Wm −2 K −4 ) ε is emissivity of the surface (-) Thermal emissivity ranges between 0 and 1 and represents the object's ability to emit radiation. A perfect emitter (i.e. black body) has an emissivity coefficient of 1 (Hillel 1998;Vollmer and Möllmann 2018), while grey bodies (ε < 1) emit less energy by radiation than a black body at the same temperature. For example, metals are usually characterized by very low emissivity (ε < 0.1), while rocks have generally high emissivity (ε > 0.9; Fiorucci et al. 2018;Mineo and Pappalardo 2021).
Infrared thermal cameras measure the infrared radiation emitted by bodies in the spectral band 7.5-14 μm. Nevertheless, in both indoor and outdoor conditions the signal received by the sensor includes contributions from the (4) J = T 4 Fig. 3 Cooling experiment setup. a Indoor setup, including a thermohygrometer to monitor air conditions (1), a rough gold panel to quantify reflected thermal radiation (2), a neoprene background (3), and a FLIR T1020 high-resolution IRT camera to monitor sample temperature (4). b Example of raw thermal image acquired during cooling experiments, that were repeated looking each sample from four perpendicular viewpoints (c). Outdoor experiments were performed on gneiss samples using a similar setup (without neoprene background) under natural solar forcing. Shaded area refers to the cooling stage (d) surrounding environment (i.e. "reflected radiation"), also modulated by the thermal transmittance of the atmosphere between the emitting object and the sensor. These components must be removed to properly interpret thermal images, or thermograms (Shannon et al. 2005). When multiple images are collected in time, time histories of temperature can be extracted from multi-temporal stacks of thermograms. IRT technique is therefore optimal to acquire thermal response data from fractured rock samples, preserving their state and physical properties making possible even high repeatability of measurements.

Cooling Experiments
We performed two series of cooling experiments as follows: (a) laboratory experiments in controlled environmental conditions under forced heating; (b) an additional outdoor experiment in natural conditions, where samples were heated by solar radiation in a diurnal cycle.
For both laboratory and outdoor experiments, we monitored the temperature of samples over time using a FLIR T1020 high-resolution thermal camera (1024 × 768 pixels), characterized by a measurement sensitivity within 0.02 °C and a field of view (FOV) of 28° × 21° (Fig. 3a). For all the experiments, the camera was mounted on a tripod located 1-2 m from the samples.
To perform cooling experiments in controlled conditions, we prepared an original laboratory setup (Fig. 3a) hosting up to three well-spaced samples for each experiment, with a high-emissivity neoprene background (ε ~ 1) to avoid disturbances around the sample due to reflected thermal radiation components. For each experiment, ambient temperature (i.e. the lower temperature boundary condition for cooling) and humidity were kept constant (20 °C and 33%, respectively) and monitored by a thermo-hygro data logger. A rough gold panel was placed in the scene to allow the measurement of reflected thermal radiation components (Fig. 3a). All the gneiss (G) and mica schist (M) samples were heated in a laboratory oven until reaching an homogeneous temperature of 80°, then placed on the set. The heating temperatures of the samples were chosen based on previous experiences on the subject Guerin et al. 2019), to embrace the temperature ranges usually encountered in natural field conditions and emphasize the cooling rates of different samples.
Sample temperature during cooling was monitored by capturing thermograms every 2 min in time-lapse mode (Fig. 3b); each indoor tests lasted approximately 1 h and 45 min. For each sample, experiments were repeated with the same initial and boundary conditions to monitor the cooling of samples from four perpendicular directions (Fig. 3c). This allowed covering the entire sample surface and capturing the effects of fracture heterogeneity on the bulk thermal behavior of studied fractured rocks.
For outdoor experiments we used a similar setup without the neoprene background, to allow samples heating under natural solar forcing over a day (Fig. 3d). Based on the results of the cooling experiments performed in the laboratory, for outdoor experiments we decided to focus on gneiss samples (G) that showed the clearest signals. Samples were arranged in a N-S trending row, exposed to solar radiation from the early morning until sunset of a cloud-free May day and then let cool for 14 h (Fig. 3d). In this experiment, thermograms were captured every 15 min in time-lapse.

Thermal Image Processing
We imported the acquired thermal images into the software FLIR ResearchIR Max™ and set all the relevant corrections constant for the entire duration of each experiment (i.e. sample-camera distance, air temperature and humidity, reflected temperature, sample rock emissivity). Both the studied rock types, that represent similar mixtures of similar mineral species (chiefly quartz, mica and feldspars) are characterized by high broadband emissivity (ε > 0.95) with values scattered over similar ranges (Salisbury and D'Aria 1992;Rivard et al. 1995;Mineo and Pappalardo 2021).
We first inspected the thermal images (Fig. 4) to check the quality of the acquired data (absence of ambient noise/disturbances, border effects related to the sample shape, etc.). Visual inspection also provided key evidence to support the quantitative interpretation of temperature distributions and temperature time histories, accounting for the influence of individual fractures (Fig. 4), their aperture (Fig. 4) and spatial distribution on the overall cooling behaviour of each sample.
For each sample, we extracted populations of temperature values from the images of the four considered sectors (Fig. 3c) using Regions-of-Interest (ROI). ROIs were selected to cover the sample surfaces excluding the external borders, where the signal was affected by edge curvature effects. We analyzed the statistical distributions of temperature in space and time using an original MATLAB™ tool, that is also able to perform a peak analysis using a Kernel distribution fitting, that automatically identify the number of frequency peaks and their maximum separation (MPD, Maximum Peak Distance, i.e. the distance between the highest and the lowest peak) in temperature distributions.
Temperature values extracted for each sample ROI at different time were averaged to obtain cooling curves representative of the bulk sample surface. We further processed these time histories to compute cooling rates in the form of the Cooling Rate Index (CRI = ΔT/Δt; Pappalardo et al. 2016). Since computing the CRI for different time intervals may hamper a significant comparison between different samples with non-synchronized cooling histories, we computed the CRI for different temperature intervals using original MAT-LAB™ script. For example, the code CRI 6070 corresponds to the CRI (in C°/min) related to the cooling curve portion between 70 °C and 60 °C.
Finally, we parametrized the lumped capacitance equation (Eq. 3) on our experimental data, to provide a description of the cooling curves in terms of its exponent, here referred to as Curve Shape Parameter (C SP ; Eq. 5): The C SP (s -1 ) provides a synthetic description of the shape of experimental cooling curves, independent on absolute thermal boundary conditions, yet incorporating thermal properties of the rock material (k,α), the boundary convective condition of the solid-air interface (h) and the representative size of radiant element (i.e. rock blocks / fragments).

Statistical Analysis
We investigated the relationships among cooling dynamics descriptors (i.e. CRI, C SP ) and fracture abundance measures (P 32 , P 33 ) by performing statistical regression analysis of the (5) C sp = h kL c experimental data. We evaluated the suitability of different regression models (linear and non-linear) by analysing residuals to check for independence of the error terms and absence of process drifts. We applied non-parametric statistics in terms of Spearman's rank correlation coefficient, ρ (Spearman 1904;Zar 2005) to evaluate the occurrence and strength of statistical relationships among variables, independently on data normality and linearity of regression models (Swan and Sandilands 1995). We also quantified the goodness of fit in terms of adjusted coefficient of determination (R 2 ) and parametric Pearson's correlation coefficient, r.

3D Finite-Element Modeling
We investigated the physical processes underlying the observed cooling behavior of selected samples using 3D Finite-Element numerical modeling, to support the interpretation of the empirical relationships between cooling dynamics and rock fracturing state. We ran 3D, continuumbased, non-linear FEM simulations of heat transfer using the code MIDAS GTS-NX™. We generated solid models with regions with properties typical of air ( Fig. 5; Table 2), that mimic the fracture geometry, size and aperture of three The analysis domains were discretized into FEM 3D meshes consisting of hybrid finite elements with size from 4.6 mm to 2.6 mm depending on the model and related fracture pattern.
Model thermal boundary conditions were set to reproduce the laboratory cooling experiments. Heat transfer simulations were conducted using a transient model accounting for both heat diffusion (in solids) and convection (at solid-fluid interfaces, i.e. at sample boundary and within open fractures). Typical literature values of the thermal properties like thermal conductivity κ Wm −1 K −1 and specific heat capacity C p JKg −1 K −1 were adopted and kept constant for the gneiss material and the air (Table 2; Barreira and Freitas 2007;Clauser and Huenges 1995). Different values of the convection coefficient h Wm −2 K −1 were applied at the solid-air interface outside the sample and inside fractures ( Table 2).
We performed numerical simulations in different steps, namely the following: (a) steady-state heating, to make every element of the model reach the same starting temperature of 80° (Fig. 5c) reproducing the oven heating in lab; (b) transient-cooling, as a time step analysis to simulate the laboratory cooling phase down to 20 °C (room temperature). We

Fracture Patterns and Thermal Response
Thermal images, acquired at different cooling stages in laboratory experiments, clearly show the influence of fracture orientation, size, aperture and intensity, typical of the two rock types, on the overall thermal response monitored by the thermal camera (Fig. 4). In fact, the complexity of the folded fabrics resulted in complex brittle failure mechanisms and strongly heterogeneous fracture networks (Agliardi et al. 2014). Moreover, unconfined conditions allowed dilatancy to occur during the laboratory fracturing process, forming fractures up to few millimetres open and leading to the development of significant fracture porosity. In the quartz-rich gneiss samples ( Fig. 1a-g), failure across foliation or along fold axial planes involved the development of dense networks of microscopic cracks within quartz grain aggregates (Agliardi et al. 2014(Agliardi et al. , 2017, resulting in a more complex and distributed fracture damage patterns than in mica schist (Fig. 1h-n).
The complexity of rock fracture patterns is mirrored by clearly different patterns in the thermal behaviour of different rock types during cooling (Fig. 4). Gneiss samples, characterized by heterogenous networks of short fractures along and across foliation and by distributed damage along cmwide brittle shear zones, show a rather "continuous" thermal behaviour (Fig. 4a, b). For this rock type we observed the expected good linear correlation between fracture porosity, P 33 , and the volumetric fracture density, P 32 (Mauldon et al. 2001;Rogers et al. 2017), with P 33 scaling with the abundance of fractures that underwent dilatancy during brittle failure (Fig. 6a). On the other hand, mica schist samples are characterized by generally higher fracture porosity associated with spaced, persistent fractures with a wide range of aperture and show a more "discontinuous" thermal behaviour than observed in gneiss samples (Fig. 4c, d). In these samples, the linear correlation between P 32 and P 33 is weaker (Fig. 6b), suggesting a more variable fracture opening of fractures transecting entire sample volumes. This suggests that adopted fracture abundance measures provide a less effective description of the fracturing state of "blocky" fractured media.

Temperature Distributions During Cooling
Samples characterized by different fracture patterns show different evolutions of their thermal states in different stages of cooling. These are captured by looking at the statistical distributions of temperature values sampled from each pixel in the thermograms ROIs at different times.
As showed in Fig. 7, at the beginning of the cooling process (i.e. samples just removed from the oven) samples tend to exhibit a unimodal temperature distribution with variance depending on rock fabric (Fig. 7a, c). When persistent open fractures filled by air are present, they are always hotter than rock and result in additional peaks with low MPD (mean peak distance; Fig. 7b, d). As the cooling process continues, different temperature distributions can be observed. Rock samples characterized by low to moderate fracturing degree and distributed damage patterns tend to keep on behaving as "thermal continua", with unimodal temperature distributions or few peaks with very low MPD (Fig. 7a, c). On the contrary, samples with high or heterogeneous fracturing behave Fig. 6 Correlations between unbiased fracture abundance indices (P 32 and P 33 ) for the gneiss (a) and mica schist (b) sample set. Dashed lines: best fitting linear models. Shaded areas: regression prediction bands (95% probability). r Pearson correlation coefficient as "thermal discontinua", with multimodal temperature distributions characterized by several, well-separated peaks that testify to the occurrence of independent radiant elements (i.e. rock blocks) separated by persistent fractures (Fig. 7b,  d). On the long term, the overall thermal states of the samples tend to equilibrate with the surrounding environment.

Cooling Dynamics and its Physical Controls
Laboratory cooling experiments show that samples with different fracturing states are characterized by different cooling curves (Fig. 8a) that have a negative exponential form according to the Newton's law of cooling (Vollmer 2009). In particular, values of the Cooling Rate Index, computed for different temperature intervals (i.e. 70-60 °C, 70-50 °C, 70-40 °C, 70-30 °C), always show a significant linear correlation with both fracture porosity, P 33 , and volumetric fracture intensity, P 32 (Fig. 8b, c; statistics in Table 3). Such correlation, that is stronger in early cooling stages, shows that more fractured rocks cool faster Mineo and Pappalardo 2016, 2022. Linear correlations are stronger for gneiss samples than for mica schist (Table 3). These results suggest that both P 33 and P 32 can be well discriminated at the accuracy level of the infrared camera during the cooling stage. Observed relationships also emerge in uncontrolled outdoor conditions (Fig. 9). In this case, the statistical correlation between the CRI values (computed over the range 40-35 °C) is weaker than those obtained in the laboratory, due to the experimental noise related to the small temperature differences and the environmental conditions (i.e. limited solar forcing during Spring, episodic change of the environmental temperature due to cloud passage or air circulation; Fig. 3d). These factors must be carefully considered when attempting to apply our experimental approach in situ, yet experimental results show its potential.
The results of 3D FEM simulations performed for selected samples (Fig. 5), processed to compute the same CRI descriptors used in laboratory experiments, provide interesting insights into the physical factors that control the cooling of fractured media. Figure 10 clearly shows that numerical modelling results reproduce very well the experimental ones, within the prediction intervals of statistical Fig. 8 Relationships between cooling dynamics and fracturing for the studied sample population. a Experimental cooling curves for gneiss (G) and mica schists (M) samples; b, c empirical correlations between the Cooling Rate Index, quantified for various tempera-ture ranges, and the fracture abundance indices (P 33 , P 32 ) for gneiss (b) and for mica schists (c). Dashed lines: best fitting linear models (Table 3). Shaded areas: regression prediction bands (95% probability) regression. This match has been obtained using constant thermal properties for the solid rock material (Table 2), while the parameter controlling convective heat transfer at solid-air interfaces (i.e. external sample boundaries and open fracture surfaces) has been tuned by calibration (Table 2). In fact, while the realistic reconstruction of the geometry of radiant elements (i.e. sample structure) was enough to guarantee a correct simulation of diffusive heat transfer, the convective parameter needed to be adapted to the different persistence and aperture of fractures in different samples.

Synthetic Description of Cooling Dynamics: Curve Shape Parameter
As showed in Fig. 8, the values of CRI depend on the considered stage of the cooling process, as well as on the initial and ambient temperatures , 2022. To isolate the control of fracturing state on the cooling behaviour of fractured samples, we expressed their temperature histories in normalized form, according to Eq. 3 (Fig. 11a) We parametrized the lumped capacitance model on normalized experimental data by non-linear regression and obtained the values of the Curve Shape Parameter (C SP ; Eq. 5) best representing the cooling dynamics of each sample. The C SP provides a unique, scalar descriptor of the shape of cooling curves and shows a good linear correlation with the unbiased fracture abundance measures ( Fig. 11b; Table 4). Also in this case, the correlation is weaker when samples with a strongly discontinuum Shaded areas: regression prediction bands (95% probability) for experimental data (Fig. 8b, Table 3) thermal behaviour are included in the analysis, as in the case of mica schist samples. In general, the C SP shows a consistent positive linear correlation with the degree of fracturing (Figs. 11b and 12a). For both the gneiss and mica schist sample sets, we extrapolated the absolute values of C SP for intact rock (C SP(int) ), i.e. the intercept of regression curve; Fig. 11b) and used them to normalize the C SP vs P 33 and C SP vs P 32 regression curves (Fig. 12b). This representation of the relationships between C SP and the fracture abundance measures seems independent on rock type and solely related to the fracturing state.

Discussion and Conclusions
A reliable assessment of the fracturing state of rocks and rock masses is key to quantify their hydro-mechanical properties and to track the processes that control the evolution of brittle rock damage at crustal levels, in natural and engineering contexts. However, this remains a difficult task, due to the intrinsically scale-dependent and heterogeneous nature of fracture networks (Dershowitz and Herda 1992;Neuman 2005;De Dreuzy 2012;Watkins et al. 2015) and to a variety of statistical biases (e.g. orientation, size) that affect fracture sampling (Priest and Hudson 1976;Baecher 1983;Mauldon et al. 2001).
The statistical sampling of discontinuities, based on field or remote mapping techniques (Priest and Hudson 1976;ISRM 1978;Sturzenegger and Stead 2009), remains a fundamental approach to the characterization of fracture networks, yet its results are affected to great extent by scale, orientation and sampling biases (Rohrbaugh et al. 2002). Non-destructive imaging techniques like X-ray computed tomography allow an effective characterization of the geometry and abundance of fractures in a small-scale three-dimensional space (Keller 1998;Agliardi et al. 2014Agliardi et al. , 2017, but have  (Table 4). Shaded areas: regression prediction bands (95% probability); r Pearson correlation coefficient. Red arrow: C SP intercept values (i.e. intact rock) (color figure online) Table 4 Goodness-of-fit statistics for the correlations between C SP and rock fracture porosity (Fig. 11 (Barton 2007;King 2009;Dorn et al. 2012;Agliardi et al. 2016) usually lack the resolution required for a detailed characterization of rock fracturing state. Infrared Thermography (IRT) applications started appearing in the fields of geoscience and geoengineering in the last 15 years, but the quantitative approaches to predict the physical properties of rocks by characterizing their thermal behaviour are few and mostly limited to the assessment of density and porosity , 2022. Our experimental investigation showed quantitatively that the cooling behaviour of low-porosity (e.g. crystalline) rocks consistently depends on rock composition and fracturing state (i.e. fracture intensity, fracture porosity and their heterogeneity) and that modern thermal cameras offer both spatial resolution and measurement accuracy suitable to detect the different thermal responses of samples characterized by different degree of fracturing.
Our experimental investigation showed that the ambivalent "continuum vs discontinuum" nature of fractured rocks, that represents the major challenge when attempting a sound characterization of their hydro-mechanical behaviour (Hoek and Brown 1997;Hoek and Diederichs 2006;Rutqvist 2015), also applies to their thermal behaviour. Recognizing thermal "continua" and "discontinua" depending on the scale of observation is the first major step towards a correct interpretation of IRT experimental data and toward a sound characterization of their degree of fracturing and fracture heterogeneity (Figs. 4 and 7). In fact, fracture abundance measures (Dershowitz and Herda 1992) like volumetric fracture intensity (P 32 ) and fracture porosity (P 33 ) provide different information for practical use, depending on whether we consider "continuum" fractured rocks, with distributed damage, homogeneous strength, deformability and permeability, or "discontinuum" materials, whose hydro-mechanical behaviour is controlled by individual, persistent fractures.

Fig. 12
Synthetic relationships between cooling dynamics and fracturing. a conceptual explanation, in a predictive perspective; b correlation between the normalized C SP and the fracture abundance indices (P 33 , P 32 ) for entire studies sample population (i.e. both G and M). Dashed lines: best fitting linear model. Shaded areas: regression prediction bands (95% probability); r Pearson correlation coefficient Our results confirmed that simple, linear correlations hold between different descriptors of cooling dynamics, as the already suggested Cooling Rate Index , 2022, and fracture abundance measures. These correlations can be established across different stages of cooling and consistently apply to sample populations characterized by different composition and fracture heterogeneity (Fig. 8). Outdoor cooling experiments suggest that observed relationship between cooling dynamics and fracturing state can be detected by IRT also in natural conditions (i.e. solar forcing and free air circulation), despite the higher experimental noise related to small temperature differences, weather and atmospheric effects.
Our 3D finite-element simulations of heat transfer, reproducing the observed cooling behaviors, allowed to support the experimental results and investigate the underlying physical processes and controls (Fig. 10). Model results suggest that, consistently with the theory, the cooling process in a fractured medium is mainly controlled by the geometry and typical size of individual radiant elements (i.e. rock blocks forming the samples), that controls heat diffusion and by the aperture of major (persistent) fractures that influence the convective component of heat transfer. According to the models, the latter has a major influence on the bulk thermal response of fractured samples and may be the responsible of the higher temperatures of open fractures in early cooling stages.
Since the observed experimental cooling behaviour strongly depends on the absolute thermal boundary conditions (initial and ambient temperatures) and by rock composition (lithology), it is generally difficult to use empirical descriptor like the Cooling Rate Index to describe rock fracturing in a practical predictive perspective, although it shows a clear correlation with rock porosity or geomechanical quality related to fracturing by previous works Mineo and Pappalardo 2016, 2022. In this work, we overcome this limitation by parametrizing a theoretical description like the lumped capacitance model (Vollmer 2009;Bergman et al.2011) on using experimental data, thus obtaining a synthetic descriptor of the shape of cooling curves, i.e. the Curve Shape Parameter (C SP ). The C SP , intrinsically referred to a continuous, homogeneous description of the fractured medium at the scale of consideration, is linearly correlated with unbiased fracture abundance measures (Fig. 11) and independent on rock type (Fig. 12). Such linear correlation, that requires some experimental quantification of fracture abundance to be established, is statistically strong (and thus suitable for prediction) when the C SP values are normalized with respect to reference "intact-rock" ones.
Our results, although not directly applicable to in situ fractured rock masses due to different fracture network characteristics on sample and outcrop scales, lay the foundations to upscale our methodology to field conditions. This would require accounting for the radiative characteristics of natural environments, the limitations of the technique and the scale effects typical of fractured rock-mass, to achieve a contactless evaluation of fractured rock mass conditions for practical applications, e.g. block size distribution evaluation, rock mass quality designation, hydro-mechanical characterization for geological and engineering application.