Quantifying the Permeability Enhancement from Blast-Induced Microfractures in Porphyry Rocks Using a Cumulant Lattice Boltzmann Method

The permeability of rocks is important in a range of geoscientific applications, including CO2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {CO}_{2}$$\end{document} sequestration, geothermal energy extraction, and in situ mineral recovery. This work presents an investigation of the change in permeability in porphyry rock samples due to blast-induced fracturing. Two samples were analysed before and after exposure to stress waves induced by the detonation of an explosive charge. Micro-computed tomography was used to image the interior of the samples at a pixel resolution of 10.3μm\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10.3\,\mu m$$\end{document}. The images were segmented into void, matrix, and grain to help quantify the differences in the rock samples. Following this, they were binarised as void or solid and the cumulant lattice Boltzmann method (LBM) was applied to simulate the flow of fluid through the connected void space. A correction required with the use of inlet and outlet reservoirs in computational permeability assessment was also proposed. Interrogation of the steady-state flow field allowed the pre- and post-loading permeability to be extracted. Conclusions were then drawn as to the effectiveness of blasting for enhancing fluid accessibility via the generation of microfractures in the rock matrix within the vicinity of a detonated charge. This paper makes contributions in three fundamental areas relating to the numerical assessment of permeability and the enhancement of fluid accessibility in low-porosity rocks. Firstly, a correction factor was proposed to account for the reservoirs commonly imposed on digitised rock samples when investigating sample permeability through numerical methods. Secondly, it validates the benefits of the LBM in handling complex geometries that would be intractable with conventional computational fluid dynamics methods that require body-fitted meshing. This is done with a novel implementation of the cumulant LBM in the open-source TCLB code. Finally, the improvement in fluid accessibility in low-permeability rock samples was shown through the assessment of multiple regions within two blasted samples. It was found that the blast-induced loading can generate extended microfractures that results in multiple orders of magnitude of permeability enhancement if the target rock possesses existing weaknesses and/or mineralisation. A novel implementation of the cumulant lattice Boltzmann method is validated for flows in complex geometries without the need for body-fitted meshing, providing an advanced tool for quantifying permeability of sub-surface samples; Quantifies the potential enhancement in fluid accessibility achieved through blast-induced damage in a confined environment through a combined workflow of micro-CT imaging, blasting experiments and high-resolution fluid simulations; Derives a correction factor to account for in-flow and out-flow reservoirs commonly introducted in simulations through digitised porous samples. A novel implementation of the cumulant lattice Boltzmann method is validated for flows in complex geometries without the need for body-fitted meshing, providing an advanced tool for quantifying permeability of sub-surface samples; Quantifies the potential enhancement in fluid accessibility achieved through blast-induced damage in a confined environment through a combined workflow of micro-CT imaging, blasting experiments and high-resolution fluid simulations; Derives a correction factor to account for in-flow and out-flow reservoirs commonly introducted in simulations through digitised porous samples.


Introduction
The permeability of a porous and/or fractured rock mass governs the percolation of fluid within a subsurface environment. It is an important hydrogeological transport parameter in a number of fields, including CO 2 sequestration (Heinze et al. 2015), geothermal energy extraction (McClure and Horne 2014), oil and gas recovery, and acid leaching of valuable minerals either in situ or within dump-or heap-leach operations (Sinclair and Thompson 2015). In all of these examples, the permeability depends on the connectivity of the void space observed in the form of pores and fractures. Due to the complex microstructural morphology of a typically heterogeneous rock mass, the permeability parameter can be timeconsuming and difficult to accurately evaluate (Holmes et al. 2015;Vianna et al. 2020).
The permeability of a rock mass is subject to change through tectonic, and other naturally occurring phenomena, but also from anthropogenic influence. Often this human influence is purposefully imposed with the aim of enhancing the flow of fluids through the porous and/or fractured environment. This is particularly true in the examples previously introduced. Hydraulic fracturing is arguably the most common form of permeability enhancement, with widespread use in the oil and gas industry (Ren et al. 2018). However, this method often propagates a single or few dominant fractures from a perforation point, which is not favourable for all applications. For example, if one is looking to leach minerals in situ, a distributed network of fractures becomes desirable to increase the fluid accessibility to the target mineral. Alternative techniques such as blasting have been discussed in the literature since the early 1980s (Adams et al. 1981). Blasting represents an efficient means of rock fracturing and fragmentation that is commonly used in conventional mining applications, and its application to generating fractures has also been acknowledged (Gharehdash et al. 2021). Further, the assessment of the permeability of blasted samples was attempted in the recent work of Gharehdash et al. (2020), which applied smooth particle hydrodynamics to simulate the impact of a blast in a Barre granite rock. A realisation of the fracture network was then generated, and OpenFOAM was used to simulate the flow field and ultimately determine the permeability of the damaged sample. On a larger scale, researchers are also interested in how the permeability of a fractured rock can vary based on in situ stresses (Latham et al. 2013;Lei et al. 2017). In a study by Latham et al. (2013), it was found that stresses can alter permeability through the propagation of cracks that can both create new flow pathways and connect existing networks. The impact of this can lead to enhanced connectivity and flow, but also create highly heterogeneous stress distributions and variable aperture through a fracture network.
In addition to blasting, fractures are also induced by the in situ stress state of the subsurface. Quantifying these fractures and the potential flow through them is of critical importance, for example, in high-level radioactive waste disposal (Chen et al. 2014). It is well understood that fractures in the subsurface from thermal and mechanical cracking can lead to the variation in hydraulic properties. Various workflows have been developed to predict this variation (Zhou et al. 2020). A fundamental understanding of the fluid transport properties in these damaged rock environments depends on the ability to capture the geometric complexity of the pore space (Sun and Wong 2018). To this end, the extraction in situ rock and characterisation through 3D imaging techniques have become commonplace for permeability investigations in the field. However, this method often requires tuning parameters associated with the pore geometry in order to align the de-stressed rock to that which was initially in situ.
Given the comparatively high cost of experimentally determining permeability, numerical methods have been employed as an alternative at varying levels of complexity (Keehm 2004). The use of numerical techniques allows for both evaluation of hydrological parameters, as well as observation and study of complex transport phenomena that would otherwise be intractable in a physical experiment. Early models often looked to predict permeability through a relation (e.g. Kozeny-Carman) to the porosity of the medium, and/ or making relations to electrical analogies (Schwartz et al. 1993). With the introduction of synchrotron X-ray microtomography (also termed micro-computed tomography, and herein referred to as µCT), researchers took to describing the pore network statistically in order to generate synthetic media that possess similar transport properties to the target rock (Hazlett 1997;Yeong and Torquato 1998). Further to this, the use of µCT also provided binary, voxelised images that could be used for direct simulation of the flow through the scanned rock sample. Numerous numerical techniques have been employed on this front including the finite element method (Ghaddar 1995;Vianna et al. 2020), pore network models (Blunt 2001;Blunt et al. 2002), smooth particle hydrodynamics (SPH) (Holmes et al. 2015;Zhu et al. 1999), and the lattice Boltzmann method (LBM) (Belov et al. 2004;Pan et al. 2006;Ramstad et al. 2010;Regulski et al. 2015;Sun et al. 2017).
In this work, the propensity for blast-induced permeability enhancement of a rock matrix is investigated through a combination of experimental blasting, X-ray tomography, and numerical flow field evaluation. This is particularly motivated by the leaching of valuable minerals in which lixiviant penetration into rock particles is critical for effective recovery. In order to determine the impact of a permeability enhancement technique, an effective workflow to quantify the flow in a sample before and after treatment was essential and is documented in this study. The remainder of this paper is structured to, first, introduce the applied methodology. This outlines the µCT imaging of the rock samples, the simulation process, and how the permeability can be extracted from simulations. In this, a correction factor is proposed to account for simulations that make use of a periodic domain with inlet and outlet fluid reservoirs. Following this, the cumulant LBM implementation and the correction factor are verified through comparison with the analytical permeability of a smooth fracture. After comparing to this analytic result, the methodology is validated with the digital rock physics benchmark presented by Andrä et al. (2013aAndrä et al. ( , 2013b. With this completed, the rock images are used to characterise the two samples before analysing their hydrological properties. Finally, conclusions are made as to the observed microfracturing and potential permeability enhancement of the rock matrix that can be obtained through blasting.

Methodology
In this study, the LBM was used to determine the permeability of various regions within two porphyry rock samples both before and after they were subjected to loading from an adjacent blast. This was achieved by first embedding a rock sample into a cylindrical grout medium with two blast holes drilled on either side of the sample to concentrate a stress wave. The samples were imaged using µCT prior to embedment and after dislodgement due to the induced blast wave. The images were then segmented and the open-source lattice Boltzmann framework, TCLB (Łaniewski-Wołłk and Rokicki 2016), was used for the simulation of fluid flow through the microfractures in the samples. The final stage of the study was post-processing to determine the permeability in each Cartesian co-ordinate to assess the impact of the blast wave. It is highlighted here that it is common practice to use a single-or multiple-relaxation time kernel for simplicity of the LBM scheme; however, the advanced cumulant relaxation scheme was employed for this investigation (Geier et al. 2015). The methodology for each of these components is discussed in this section.

Experimental Blasting Design
The two grout cylinders had a diameter of 230 mm and a height of 60 mm . Both the samples had two blast holes drilled with diameters of 8 mm . The blast hole locations are indicated in Fig. 1 by the distance d and contained 0.725 g of PETN explosive. For Sample 1, this distance was 44 mm and for Sample 2 it was 54 mm . The dimensions of the two samples are supplied in Table 1. Figure 1 provides a schematic of the blasting experiment in order to generate microfracturing in the embedded samples. The medium used to surround the rock sample was a grout acquired from Fosroc known as Conbextra GP. 1

Imaging of Rock Samples
Two semi-cubical porphyry rock samples of different size, both before and after being subjected to stress induced by a blast (see Fig. 2), were analysed by µCT. The samples were positioned, and stabilised with foam padding, in 34 mm tubes before being scanned with an isotropic voxel size of 10.3 m 3 . The samples were scanned in air at an energy of 90 kVp and a current of 155 A . The integration time was set at 1320 ms , once averaged, resulting in a 1.32 s sample time. A 0.1 mm copper filter was used. The scan settings were kept constant for all the analysed samples. The scanning parameters were also maintained at similar values. The images were exported as a standard Digital Imaging and Communications in Medicine (DICOM) stack to allow for further processing. The resulting µCT dataset   contained four sets of scans with resolutions of 10.3 m per pixel. The size characteristics of the µCT dataset are listed in Table 2. The µCT images were segmented using a combination of watershed segmentation (Amankwah and Aldrich 2010), active contouring (Hemalatha et al. 2018) and Sobel filter maps (Yao 2016). Figure 3 shows the three main classes that were determined, namely, grains, voids and matrix. Grains were segmented using the watershed method, and the Sobel filter was used to determine their boundaries. Voids, including pores and fractures, were extracted by the watershed algorithm. Active contouring was implemented to detect the edges of voids. Two high-resolution sets ( 2 m per pixel), obtained for smaller Examples of µCT data before and after segmentation from the Sample #1 high-resolution set, showing a the 3D volume before segmentation, b the 3D volume after segmentation, c a 2D image before segmentation and d a 2D image after segmentation. On the images before segmentation, white corresponds to grains, grey is matrix, and black represents pores. On images after segmentation, red corresponds to grains, grey is matrix, and blue represents pores. The grains and matrix are treated as solid in the LBM simulations, but provide scope for reactive transport modelling in future studies sub-volumes of the same rock samples, were used to calibrate the results of segmentation for the presented scans. Additionally, the segmentation of high-resolution images was validated by comparison with scanning electron microscopy (SEM), which was performed on the sub-volumes cut from the rock samples. In order to conduct SEM, each sub-volume was attached to a solid cylinder stub ( 25 mm in diameter) and scanned with electron microscope. The equipment used a Schottky (Thermionic) field-emission (T-FE) electron gun.
The T-FE electron gun can generate stable probe currents from pico-amperes (pA) through to one hundred nano-amperes (nA). Each sub-volume was analysed from three orthogonal sides in secondary and backscatter electron modes with an accelerating voltage of 20 keV . The working distance was set to 10 mm , and spot size of eight was chosen. SEM images have a pixel size (nominal resolution) in the range between 0.05 and 3.8 m . SEM images were segmented using the Otsu algorithm (Otsu 1979). It is highlighted here that SEM only provides imagery of the surface of a sample, and as such cannot be used to create a 3D reconstruction. In this study, the surface image obtained by SEM was used to validate the segmentation parameters for the relevant surface of the µCT scan before these were applied to the entire 3D volume.

Governing Equations and Numerical Solution Methodology
In order to resolve the flow through the domains constructed in this study, the incompressible Navier-Stokes equations, are resolved using a cumulant LBM. Here, is the macroscopic velocity at a given location and point in time, p is the pressure, is the kinematic viscosity, and is the fluid density. As previously stated, the LBM was used to recover these hydrodynamic equations.
The LBM was used for the solution of the Navier-Stokes equations. A good introduction to the LBM can be found in Succi (2019) and Krüger et al. (2017). The method makes use of a Cartesian lattice of equal, non-dimensional spacing, Δx = 1 , and a fixed non-dimensional time step, Δt = 1 . On the defined lattice a set of streaming velocities, { e 1 , … , e q }, is selected to match the grid spacing. For the discretisation, a set of density fields, { f 1 , … , f q }, is given in each point of the lattice, from which macroscopic quantities can be determined as, where c s is the speed of sound. In this study, the so-called D3Q27 lattice is used, which is a three-dimensional Cartesian lattice, with velocities, e, spanning all 26 first-order neighbours, and the zero-velocity vector. The D3Q27 lattice, along with the D3Q15 and D3Q19, is the most common discretisations of the velocity space to resolve the Navier-Stokes equations in LBM literature. Although the lattice used here requires more memory and computation with streaming directions to all neighbour nodes, it provides higher accuracy (Krüger et al. 2017). The speed of sound on this lattice is equal 1∕ √ 3. The iteration of the numerical scheme consists of two components, namely, streaming and collision. The streaming step moves the densities between the nodes of the lattice, while the collision step describes the application of a nonlinear collision operator, W , to densities in all nodes of the lattice. The combined iteration can be expressed as, One of the most commonly used collision operators is the Bhatnagar-Gross-Krook (BGK) (Succi 2019), , and f eq is the equilibrium distribution. The equilibrium distribution function is found by matching the moments of the distribution to the moments of the continuous Maxwell-Boltzmann equilibrium (normal) distribution (Succi 2019), Alongside the equilibrium distribution, most of the recent advancements in the LBM can be conveniently expressed in moment space. The raw moments of the distribution can be defined as, For an appropriate selection of q moments, one can define an invertible transformation, m = Mf , which computes this vector of moments from the vector of densities. In turn the collision operator can be expressed using moments, which is frequently more computationally efficient, as moments of equilibrium are faster to compute then the equilibrium distribution itself. In addition to this, one can choose different relaxation factors, , for certain moments, resulting in what is called the two-relaxation time (TRT) (Ginzburg 2005), or multiple-relaxation time (MRT) collision operator (d 'Humières 2002).
The present paper uses the cumulant LBM, which in place of moments uses cumulants. As cumulants are a lesser-known concept of the two, a brief description is included in this work. For a more detailed discussion of the method, the reader is referred to the work of Geier et al. (2015). For any random variable, X, one can expand the natural logarithm of the momentgenerating function into a Taylor series, Cumulants of X are defined as the coefficients of this expansion, c i . It is easy to show that c 0 , c 1 and c 2 , are equal to zero, the mean and the variance of X, respectively. This definition can be easily extended to n-dimensional random variables. Moreover one can define an extension of the definition to distributions, which do not integrate to 1, by normalisation. As there are well-known equations to compute cumulants from raw moments, and vice versa, one can define a transformation, c = Cm , which computes the vector of cumulants from the vector of moments. This vector of cumulants has to additionally include the zeroth moment, namely the sum of the f, for the operator to be invertible. This allows one to express the collision operator using cumulants as (Geier et al. 2015), where cumulants of the equilibrium distribution, c eq , are more straightforward to compute than the equivalent moments, m eq . This is because any normal distribution has third-and higher-order cumulants equal to zero. Again, different relaxation factors can be used for different order cumulants in order to enhance the stability of the collision operator. It is important to note that in contrast to M , the cumulant transformation is not linear; thus, even for a single choice of relaxation factor, , for all order cumulants, it is not equivalent to the BGK collision operator.
In this work, all cumulants are relaxed with = 1 , apart from the second-order cumulants, which are relaxed based on the kinematic fluid viscosity, = 1 3 + 1 2 . Moreover, the sum of cumulants, c 200 + c 020 + c 002 , is relaxed with = 1 . More details can be found in (Geier et al. 2015). The source code for the specific cumulant collision operator is openly available. 2 The relaxation factors are chosen to ensure stability, while maintaining the required degree of accuracy (Geier et al. 2015).

Characterisation of the Flow
When assessing the flow of fluid through a porous medium, the absolute permeability, k, is a key parameter of interest (Holmes et al. 2015). This metric provides a quantitative measure of how easily a fluid will permeate through the porous medium of interest. In this work, the absolute permeability will be specified in milli-Darcies (mD).
In this study, the method described by Holmes et al. (2015) was employed to determine the absolute permeability of test samples. In order to do this, one first starts from Darcy's law given as, where is the dynamic viscosity of the saturating fluid, U is the superficial fluid velocity, and ∇p is the mean driving pressure gradient. To solve for the absolute permeability, both U and ∇p must be determined from the numerical simulation of flow through the target sample.
Numerical simulations allow for the interrogation of local flow features, allowing for the mean flow velocity, u, to be integrated at a desired time as, where ( ) is the local density at location , u i is the local fluid velocity in the direction i, and Ω represents the flow domain, which can be a slice of the full domain due to continuity (Holmes et al. 2015). The superficial velocity can then be recovered as, in which = V f luid ∕V total is the porosity of the domain. According to Larson and Higdon (1989) the mean pressure gradient for flow over an inclusion is simply given by the force on that inclusion, divided by the volume of the unit cell in which it is contained. Here, the drag force, F d , is determined as a balance of the momentum before and after the relaxation process in the solid boundary cells, leading to, In this study the fluid phase was driven with a gravitational acceleration, . As such the pressure gradient should be equal to , providing a method to check the accuracy of the drag force calculation. The calculated permeability can then be given by, In practice, care needs to be taken with the drag force determined from a simulation. This becomes evident as it is commonplace for inlet and outlet reservoirs to be applied to ensure periodicity and reduce flow disturbance. As such, if one naively determines the drag force felt by all the solid surfaces in the domain, the additional hydrostatic pressure caused by this region will reduce the accuracy of the predicted sample permeability. Here, two methods to account for this error are introduced. The first method is to simply extract the drag force on the solid surfaces within the sample of interest, neglecting those bounding the inlet and outlet reservoirs. However, if this is not possible (e.g. experimental or software restrictions), a correction factor can be derived to reduce the associated prediction error and remove dependence on the reservoir length.
This correction assumes that the pressure drop across the reservoir is equal to zero. Thus, if a constant pressure gradient, ∇p , is applied across the domain, the gradient acting on the sample, ∇p samp. , can be calculated by, where L is the length of the sample, and l is the length of the reservoir. This gives an equation for the corrected permeability, The assumption of zero pressure gradient over the reservoir is a simplification, as the pressure drop across an outlet to an open reservoir is nonzero, and is velocity dependent. Nevertheless, as is shown in Sect. 3.1, this simple correction, which will be called the Łaniewski correction factor in this paper, can give a substantial improvement in accuracy and consistency.

Permeability Assessment of a Simplified 2D Channel
As discussed in Sect. 2.4, the simulation of porous media flow through scanned rock segments often uses a reservoir section to eliminate the issue of aligning fluid cells for periodicity. Without such a reservoir, the rock sample would need to be mirrored to ensure flow path continuity, but the computational cost of simulating twice the desired domain can be infeasible for larger test cases. In LBM, it is also standard practise to impose a body force on the domain to create a pressure gradient. This section looks to quantify the effect of the reservoir by comparing the additional pressure gradient contribution of the fluid in the reservoir sections themselves and provides practitioners with two methods to mitigate its effect. This was conducted through assessing the permeability of a fracture with aperture, h, that is calculated as, k f = h 2 ∕12.
The flow domain was constructed as a pseudo-2D Hele-Shaw cell with a length of L = 800 m , a height of H = 20 m , and a reservoir with four times H and a length of l ∈ {0, 40, 80, 120} m . The simulation resolution was constructed such that H was discretised at four different resolutions, namely {16, 32, 64, 128} cells. The simulations were periodic in all directions with the depth set to four cells. A 2D schematic of the domain is provided in Fig. 4.
In order to define the fluid properties, an approximate Reynolds number was defined and specified at 0.01 to be sure there would be minimal impact from flow separation and inertial influences at the inlet and outlet of the channel section, Here, it can also be seen that the kinematic viscosity has been set to 10 −6 m 2 ∕s to represent a fluid similar to water. Using the expected velocity from Poiseuille flow and neglecting the inlet and outlet reservoirs allows the average velocity to be related to the required gravitational acceleration, g, to achieve the desired Reynolds number, A summary of the physical and lattice parameters is provided in Table 3 for the interested reader.   Figure 5 presents an example of the steady-state total velocity in which a high velocity is observed in the channel compared with the slow moving reservoir regions. The fluid can be seen to expand as it exits the channel before converging back to the inlet region. The maximum observed velocity for this case is approximately 0.56 mm∕s , which is consistent with the initial predictions used to specify the fluid parameters for the simulation. Figure 6 shows the error associated with the prediction of the channel permeability when considering only the internal surfaces of the flow domain. Here, it is clear that the simulated permeation of the fluid through the channel and the procedure to quantify the permeability is accurate. The mean error neglecting the reservoir size was 0.28%, 0.20%, 0.18% and 0.17% for resolutions of the channel height of 16, 32, 64 and 128, respectively.
If the external solid surfaces are included when calculating the drag force in the simulation, the error significantly increases. Figure 7 presents the error in the fracture permeability prediction as the resolution of the channel increases for various fluid reservoir lengths without using the Łaniewski correction factor. Here, it can be observed that the solutions are converging towards a permeability higher than the expected value, as the hydrostatic  Error associated with the fracture permeability prediction for varying fracture resolution and inletoutlet reservoir length when calculating the drag force from interior solid surfaces pressure gradient in the reservoir is not accounted for. As such, the channel permeability prediction becomes highly dependent on the inlet and outlet reservoirs. However, this is the permeability of the entire simulated domain, and not the test sample itself. The Łaniewski correction factor in Equation (18) was derived to remove this dependence. Figure 8 presents the results of permeability predictions for all cases taking into account the Łaniewski correction factor. Here, it can be observed that each case converges more consistently when compared with the results of simulations neglecting the correction term. It is noted that the level of error is still significant for low-resolution channels; however, it is the consistency of behaviour that the correction factor aims to improve.

Permeability Prediction of a Fontainebleau Sandstone
In this section, the numerical methodology was validated by applying it to the Fontainebleau sandstone benchmark case discussed by Andrä et al. (2013aAndrä et al. ( , 2013b. The sandstone sample image was initially acquired by Exxon-Mobil using a synchrotron source at Brookhaven National Laboratory. The data were provided in the supplementary material of Andrä et al. (2013a) as a segmented binary image ( 288 × 288 × 300 voxels) with each voxel characterised as a grain or a pore. For the sample, the voxel edge length was given as 7.5 m. In Andrä et al. (2013b) the LBM and Explicit Jump (EJ) method were applied to simulate flow through the rock sample to assess the single phase, absolute permeability. It was reported that the EJ method consistently over-predicted the result of the LBM simulation, and the Fontainebleau sandstone test was no exception to this. The interested reader is directed to the works of Andrä et al. (2013aAndrä et al. ( , 2013b for further details on the methods applied by these researchers. Figure 9 presents a view of the Fontainebleau sandstone dataset with a binary segmentation as well as the reference axis used in this work. It is noted that the x-axis and y-axis are switched in comparison with the supplementary material results supplied by Andrä et al. (2013b), simply due to the method used to import the binary image into the TCLB flow code (Łaniewski-Wołłk and Rokicki 2016). Gomez (2009) produced permeability versus porosity data for a Fontainebleau sandstone that agreed well with data available in the literature and was cited as a validation source by Andrä et al. (2013b). The image displayed in Fig. 9 was extracted from a larger sample which was reported to have a porosity of 15.2 % and a permeability of approximate 1100 mD . It is noted here that the numerically tested sample in this work is a subset of the larger sample, and as such, this value simply gives an order of magnitude for the permeability of the sub-volume used. The results shown here are not aimed at indicating that the cumulant LBM is more accurate than existing methods. In this section, the aim is to show that the relaxation operator maintains the accuracy of Fig. 9 The Fontainebleau sandstone dataset in which red (online) represents solid and blue (online) represents void space. The total image size is 288 × 288 × 300 voxels, each with a side length of 7.5 m commonly used methods while gaining its associated benefits. These benefits are primarily associated with a wider range of parameters for which the method will produce stable results due to, for example, its stronger dissipation of spurious pressure waves (Geier et al. 2015).
In the numerical simulation of this sample, the cumulant LBM introduced in Sect. 2.3 was used to model fluid propagating through the void space. Inlet and outlet reservoirs were applied in order to assist with periodicity, each having a length of 16 lattice cells for all flow cases (i.e. a total reservoir length of 32 lattice cells). Flow was driven through the sandstone sample with a body force applied to the full domain in the direction of the three primary axes and as such three simulations, one for each flow direction, were conducted. During the simulation, an integral of the momentum in the direction of the applied body force was calculated every 100 iterations. This was used as a convergence criterion for the simulations with a change limit of 10 −8 placed on the momentum in lattice units. Simulations were terminated if the change in momentum was less than this limit for five sequential integral calculations. With this stop criterion, the simulations were run for 259,000, 241,000 and 223,000 iterations for flow in the x-, y-and z-directions, respectively. The simulations were run on a single Nvidia V100 GPU on which they achieved approximately 1575 million lattice updates per second (MLUps) and required 12 GB of memory. This update speed is equivalent to ∼685 GB/s of data access, which is 75% of the bandwidth of this graphics card's memory (900 GB/s). As the LBM is frequently a memorylimited scheme, the memory bandwidth can be used as an indicator of the theoretical maximum performance. It is important to note that the code is designed for use on multi-GPU architectures and achieves near perfect scaling (Fakhari et al. 2017), which allows this methodology to be applied in much larger geometries. Table 4 summarises a comparison between the simulated results of Andrä et al. (2013b) and the LBM results generated in this study. It is evident that the implementation of the cumulant LBM in the open-source code, TCLB, produces results that strongly align with the benchmark study. Further to this, the result validates that the methodology to determine permeability values is comparable with expected results for real rock samples. Figure 10 presents two views of the void space in the Fontainebleau sandstone with pressure and velocity magnitude indicated by the colour. These images indicate the nonlinearity observed in the flow within porous samples, with regions of high pressure and pore throats with relatively large velocity seen throughout the domain. The digital rock physics workflow allows the flow through these regions to be resolved and accounted for when assessing the permeability of real materials.
With the verification and validation of the permeability assessment methodology completed, the remainder of this manuscript describes the analysis of a porphyry rock sample both before and after it is subjected to loading from a blast wave. These two scenarios will be referred to as before stress and after stress in the remainder of this work. Insight into the homogeneity of the permeability tensor in addition to the expected magnitude of flow enhancement provides information, for example, to in situ leaching and or cave mining operations to assess the feasibility of blast preconditioning prior to extraction.

Imaging Analysis
The µCT images before and after blasting were aligned and then segmented to identify each pixel as rock, mineral or void. From here, three regions from each rock sample were selected in order to study the permeability evolution induced from the blast wave interaction. Figure 11 presents the results of a preliminary study of the before blast sample which demonstrated that both mineral grains and voids are scattered heterogeneously; however, their dissemination follows the patterns which reflect the metasomatic processes that the rocks underwent (e.g. sericitisation). Visual examination of more than 50 samples revealed that these patterns could result from healing of pre-existing rock fractures with a hydrothermal solution as part of a metasomatic alteration of the rock mass. The detailed examination of this preliminary finding is outside the scope of this current work; however, this hypothesis was used in selecting the regions of investigation for permeability analysis. If the initial supposition was correct, direct stress application might cause the formation of a system of cracks that follow the pre-existing fractures. Damage along these existing weaknesses could significantly enhance the permeability of the sample after applied stress. Thus, the regions of investigation were chosen from the samples before stress and those regions corresponded to the areas where indications of healed fractures could be observed. Figure 12 indicates that the section chosen for investigation-based off the before stress sample coincided with the most fractured areas of the after stress sample. In addition to the preliminary assessment of void and grain distribution, a detailed analysis of their location and size was also performed. Four parameters were examined including the volume of each grain, the volume of connected void space, the distance between voids and grains, as well as the distance between grains and voids.
To estimate the volume of grains and voids, the void and the grain ROIs (regions of interest) were converted into multi-ROIs so that every separated object represented its own ROI. From this, the volume of every multi-ROI was calculated with the Dragonfly software. The results for rock sample 1 are summarised in Fig. 13. Note the logarithmic y-axis indicating the dominant majority of grains ( ∼99%) having an effective radius of less than 135 m. The pores are noticeably smaller with the majority ( ∼99%) having an effective radius of less than 90 m.
The size of grains and pores is similar in Sample #1 and Sample #2; however, the quantity of grains in Sample #2 is significantly lower. Figure 14 graphs the size distributions for the grains and pores of the second sample. Here, it is clear that the number of grains was approximately 45 % less than in Sample #1, while the number of voids increased slightly. The distance between voids and grains was estimated using the concept of distance maps. Distance maps look like grey-scale images, where each voxel outside the region of interest is replaced with a grey value equal to that voxel's distance from the nearest voxel in the ROI. Dragonfly produces minimum, maximum and mean intensity values that indicate the minimum, maximum and mean distance found for each object in the multi-ROI, respectively. It could be interpreted in terms of the distance between voids and grains, as shown in Fig. 15. To estimate the distance between voids and grains and between grains and voids, the minimum intensity was used, which gave the distance between a void and the closest grain (for voidgrain distance) and the distance between a grain and the closest void (for grain-void distance). The results of this analysis are summarised in Figs. 16 and 17 for Sample #1 and Sample #2, respectively. Here, it can be seen that the range of separation distances is larger for the void-grain measurement. This indicates that the grains occur in more localised regions, while the pores are distributed throughout the sample. As a result of this, the grains tend to have nearby pores, while a pore can have larger separation distance if it is further from the grain regions.

Permeability Assessment of Conditioned Rock
With the use of the imaging techniques discussed in Sect. 2.2, two porphyry rocks samples were examined both before and after they were subjected to a blast wave (i.e. conditioned). Here, the focus is not on the blast experiment, but the increase in permeability within the samples. The experiment consisted of two blast holes drilled into a synthetic grout medium on either side of the embedded porphyry sample. In order to assess the permeability, three volumes of 480 and three of 768 voxels 3 were taken from rock Sample #1 and #2, respectively. This was conducted for the samples both before and after conditioning. This resulted in a total of 12 samples, six before blasting and the respective six after blasting, to be interrogated using the open-source lattice Boltzmann code, TCLB (Łaniewski-Wołłk and Rokicki 2016). The voxel length in each of the samples was 10.3 m , resulting in a cubic region of interest with a length of 4.944 mm and 7.910 mm for Sample #1 and #2, respectively. Inlet and outlet reservoirs of 16 lattice units were used for the simulations, and the methodology previously validated for permeability assessment was applied. The remainder of this section progresses through the analysis of each sub-volume looking at their state before and after loading in order to conclude the potential value of blast conditioning of rock as a method to enhance permeability.

Rock Sample #1
After imaging, Sample #1 was found to contain a lower porosity than the second sample. Table 5 lists the porosity before and after blast conditioning. Here, it is evident that the sample was significantly damaged through its interaction with the blast wave as the overall porosity increased by a factor of approximately 2.1. It is noted here that the sub-volumes were selected based on the appearance of mineral clusters as this was hypothesised to present a region of weakness in the rock where damage was expected to be induced. It is not conclusive, however, from the results of this sample, but it appears that the regions of mineralisation have behaved in this way. Often in the case of industrial applications, such as dump or heap leaching, it is not the porosity but the ability of fluid to pass through the porous medium that is important. As such, flow simulations were conducted to extract the permeability of the porphyry rock before and after blasting.

Before Blasting Permeability Analysis
Prior to launching flow simulations, the sub-volumes were first subjected to a flood fill algorithm to examine the existence of a connected flow path between the inlet and outlet regions. From this, short simulations were conducted for the sake of visualisation if no connected path was evident, or full simulations were launched if a connected path was found. Figure 18 indicates that there was no connected flow path through any sub-volume of Sample #1. Note that the images in this figure represent the 10 largest connected pore spaces to provide ease of visual inspection for connected flow paths. One could simply extract the largest connected pore space from the flood fill algorithm for this purpose; however, 10 was selected here to give the reader a view of the unconditioned sub-volumes for comparison with the conditioned analogues.
Even though there was no connected flow path observed in the sub-volumes, this does not mean the permeability of the sample is zero. The authors highlight that the Fig. 18 Flow simulation in a sub-volume 1, b sub-volume 2 and c sub-volume 3 of the before blasting rock image in which the colour represents the pressure field after 100 iterations. Simulations were terminated as no connectivity between inlet and outlet reservoir was observed at this scale scanned resolution was 10.3 m per voxel, and as such, it is expected that smaller features through which fluid could travel under either pressure or capillary action would allow for transport through such a structure. However, the lack of connectivity at the scale examined does imply a relatively low sample permeability, expected to be in the nano-Darcy range based on the results of conditioned sub-volumes.

After Blasting Permeability Analysis
Following the blast, the conditioned ore was imaged and the sub-volumes aligned in order to extract the same representative region as examined in the before blast simulations. Figure 19 provides a view of the sub-volumes in which the degree of damage is clear through the increased porosity, and development of major fractures. In particular, the sub-volumes 1 and 2 show fractures that orientate along the Cartesian co-ordinates that induce a preferential flow path when acceleration is applied in the respective direction. In comparison, sub-volume 3 shows clear damage, but the fracture curves from the y-z plane to the x-y plane. The specific Fig. 19 Flow simulation in a sub-volume 1, b sub-volume 2 and c sub-volume 3 after conditioning in which the colour (online) indicates the relative pressure in the pore space, with blue implying a lower pressure to red surface geometry of these fractures is valuable information that has been assigned to future analysis to assess the impact of tortuosity observed in realistic fractures. Table 6 summarises the permeability results extracted from the lattice Boltzmann simulations of flow through the pore space of the sub-volumes from the conditioned Sample #1. This highlights the effect of the induced fractures with significant heterogeneity observed in all sub-volumes and directional permeabilities ranging from milli-Darcy to Darcy (neglecting the z-directional flow of V2 in which no connecting flow path from inlet to outlet existed).
For this sample, it is clear that blast conditioning was able to induce significant damage in the micro-structure of the rock, enhancing porosity and local permeability by orders of magnitude. Further to this, the permeability increase appeared to come from the development of fractures through regions of high mineralisation which is a useful outcome for example to rejuvenate existing dump and heap leach operations. This conclusion was drawn based on the observed regions of mineralisation identified in the µCT images in comparison with the fractures clearly seen in Fig. 19. Future work will experimentally examine the entire sample with the use of capillary rise testing to observe this behaviour at a larger scale, as well as capture the effect of rock features below the scanned resolution.

Rock Sample #2
Sample #2 was found to have a much higher porosity prior to blasting, when compared with Sample #1. This potentially indicates a different origin or series of processes to which the rock had been exposed prior to delivery for the purposes of this study. The variation in rock provides an interesting contrast to Sample #1, in that the blast wave only generated a marginal increase in porosity. Table 7 presents the porosity for each sub-volume as well as the entire Sample #2 rock. Note here that the size of the region of interest (i.e. each sub-volume) has been increased to 768 3 pixels. In the simulation of flow through these samples, a reservoir region was again applied with a length of 32 lattice units.
Referring to Sect. 3.3, it was found that Sample #2 consisted of significantly fewer mineral grains, which appeared to be points of weakness in Sample #1. Further, it is hypothesised that the higher porosity may have acted to attenuate the blast wave at a faster rate, contributing to the reduction in damage observed for this sample.

Before Blasting Permeability Analysis
Given the higher initial porosity of Sample #2, connectivity was found in three of the nine potential flow cases (i.e. three sub-volumes and three primary flow directions for each). For sub-volume 1, a minor connectivity path was observed in the y-direction with the flood fill algorithm propagating from inlet to outlet. Additionally, connectivity in the z-direction was found for sub-volumes 2 and 3. Table 8 summarises the permeability obtained by driving flow through these domains with a body force applied in the desired direction. It is clear

Fig. 20
Flow simulation in a sub-volume 1, b sub-volume 2 and c sub-volume 3 of the before blasting rock image, in which the colour represents the pressure field after 100 iterations. Simulations were terminated as no connectivity between inlet and outlet reservoir was observed at this scale, in the x-direction of flow that although there was a connected path from inlet to outlet, the resulting permeability estimate for these is relatively low when contrasted with Sample #1. Figure 20 presents a 3D image of the sub-volumes constructed for flow in the x-direction. Here, the 10 largest connected regions are shown to indicate the topology of the pore space. The images are coloured by the calculated pressure field after 100 lattice iterations, to provide visual contrast for the flattened images. It could be argued that greater connectivity could be found with a higher-resolution scan; however, the permeability would still be expected to be in the nanodarcy range.

After Blasting Permeability Analysis
The pore space of Sample #2 after blasting can be seen in Fig. 21. Additional void and partial fractures can be observed; however, the material response is clearly different to Sample #1. This shows the impact of material heterogeneity when considering how one may effectively generate damage to the rock matrix in order to enhance its permeability. Sample #2 exhibited a higher porosity, but lower mineralisation providing a hypothesis that Fig. 21 Flow simulation in a sub-volume 1, b sub-volume 2 and c sub-volume 3 of the after blasting rock image, in which the colour represents the pressure field after drag and momentum have reached a steady state. These particular simulations were conducted with flow in the x-direction pores act to attenuate the blast wave in comparison with the mineral grains which act as a defect along which a fracture can propagate. Table 9 summarises the permeability results as determined from the cumulant LBM simulations. Here, the increased connectivity of sub-volume 1 allowed for enhanced flow in all directions. This is particularly evident in the z-direction, where an almost continuous fracture is conductive (observable in Fig. 21a). In comparison, sub-volumes 2 and 3 were found to have only minor improvements in connectivity and as a result, permeability. When examining the before blasting state of these samples, it is evident that sub-volume 1 had the lowest porosity. The increased permeability observed for this sample when compared to sub-volumes 2 and 3 supports the hypothesis of porosity-based attenuation leading to a reduced potential for fracture generation; however, further testing is proposed to examine the relevant phenomena.
Overall, this component of the study consisted of 18 simulations conducted across six flow domains and with flow driven in the three Cartesian directions. The generation of microfractures in Sample #1 led to significant enhancement in void connectivity through the sub-volumes and, ultimately, increases of multiple orders of magnitude in permeability. In contrast, Sample #2 had a higher initial porosity but saw less impact from blast wave interaction. The first sub-volume assessed for this sample was found to have a partial fracturing allowing for an increased permeability in all directions; however, sub-volumes 2 and 3 were found to have only marginal improvements.

Conclusions
This study analysed the accessibility of fluid to flow through tight, porphyry rock samples both before and after they were subjected to dynamic loading from a blast wave. Prior to this, work was first conducted to establish a correction term to allow for inlet and outlet fluid reservoirs to be accounted for in the numerical estimation of permeability. Inclusion of these reservoirs allows for periodic boundaries to be applied in the direction of flow without the challenge of inlet-outlet connectivity. However, the use of these reservoirs creates an additional pressure component that must be included when calculating the actual pressure drop across the test sample. Following this, the cumulant lattice Boltzmann implementation for permeability calculation was validated against the digital rock physics benchmark of flow in a Fontainebleau sandstone. Close agreement was found with previous numerical and experimental works for this type of porous medium.
With the correction term formulated, and the numerical method validated, detailed information was provided for the two porphyry rock samples used in this study. The first sample was observed to have a lower porosity, but higher levels of mineralisation. The results of the modelling campaign, consisting of 18 simulations, were presented in which the sub-volumes of the first rock sample indicated significant increases in permeability after being exposed to the blast wave. In contrast, the second sample showed only marginal enhancement of fluid transport. These findings led to the hypothesis that the existing pores acted to attenuate the blast wave, whereas the mineral grains provided defect points that encouraged the generation of connected micro-fractures. Further study is proposed to support this hypothesis as the outcome could be of relevance, for example, in heap, dump and in situ leaching of minerals as well as in locating potential sites for CO 2 storage. It is expected that gaining additional information on rock samples in regard to their origins and formation processes could highlight the contrasts in microfracturing between samples.