Kirill M. Gerke1,2,3,4,*, Marina V. Karsanina1,3,4,* & Dirk Mallants1

1CSIRO Land and Water, Glen Osmond, PB2, SA 5064, Australia. 2The University of Melbourne, Department of Infrastructure Engineering, Parkville, VIC, 3010, Australia. 3Institute of Geosphere Dynamics of the Russian Academy of Sciences, Leninsky prosp. 38/1, Moscow, 119334, Russia. 4Institute of Physics of the Earth of Russian Academy of Sciences, Bolshaya Gruzinskaya 10, Moscow, 107031, Russia.

Spatial data captured with sensors of different resolution would provide a maximum degree of information if the data were to be merged into a single image representing all scales. We develop a general solution for merging multiscale categorical spatial data into a single dataset using stochastic reconstructions with rescaled correlation functions. The versatility of the method is demonstrated by merging three images of shale rock representing macro, micro and nanoscale spatial information on mineral, organic matter and porosity distribution. Merging multiscale images of shale rock is pivotal to quantify more reliably petrophysical properties needed for production optimization and environmental impacts minimization. Images obtained by X-ray microtomography and scanning electron microscopy were fused into a single image with predefined resolution. The methodology is sufficiently generic for implementation of other stochastic reconstruction techniques, any number of scales, any number of material phases, and any number of images for a given scale. The methodology can be further used to assess effective properties of fused porous media images or to compress voluminous spatial datasets for efficient data storage. Practical applications are not limited to petroleum engineering or more broadly geosciences, but will also find their way in material sciences, climatology, and remote sensing.

Assembling spatial information across a wide range of scales is a crucial component in almost any type of industrial or scientific activity1. Spatial information –discrete property values with known coordinates – is usually reported by means of digital images where values are represented by intensities and coordinates are marked within a grid with known spatial sizes or image resolution (pixels for 2D images and voxels for 3D). Image resolution - a key attribute of all spatial data- has been increasing considerably in all scientific fields. For instance, performance of computed tomography has doubled approximately every two years since the mid 1980s2.

Likewise, electron microscopes (SEM) typically have a 250–300 times higher resolution than optical microscopes, while the highest spatial resolution of aerial imagery is up to 2.5 cm compared to 50 cm for satellite images. These differences in resolution by no means suggest that lower-resolution methods need to be abandoned; one of their advantages is their much wider field of view or measurement support and, thus, their ability to resolve larger objects than the more recent high-resolution devices. This becomes especially relevant for 3D imaging methods such as X-ray micro-tomography (XMT), focused ion beam–scanning electron microscopy (FIB-SEM) or neutron tomography, where the very small measurement support of the device is prohibitive of determining relevant properties across many locations of the study object. The duality between a sufficiently large, in terms of a representative elementary volume (REV), sample size and a sufficient spatial resolution can be resolved by studying the same object at different spatial scales, followed by upscaling and/or downscaling information into a single image capturing data from all scales.

Merging information across multiple scales is of interest to many scientific disciplines; some meaningful applications relate to (non-exhaustive list): 1) images of outcrops displaying large-scale features of sedimentary units, medium-scale sedimentological facies heterogeneity and fault zones, and small-scale fractures3,4 to assemble multiscale information on fracture size and connectivity distribution for simulations of flow and transport; 2) large-scale 2/3D seismic images and 1D high-resolution well-logging supported by small-scale core measurements of fluid properties to create realistic oil/gas reservoir models5 for improving resolution of reservoir characteristics and properties upscaling schemes; 3) medical imaging6 for better diagnostics; 4) merging hydrogeological parameters from aquifer pumping tests, tracer tests and undisturbed core measurements;7,8 5) fusing micrometer scale XMT and nanometer scale FIB-SEM images to derive effective properties such as apparent permeability or gas desorption rates of low-permeable carbonates, shales and coal seams9 using pore-scale modelling approaches;10 6) combining satellite images with plot-scale and point-scale on-ground measurements of plant species distributions, CO2 or CH4 emissions, or soil moisture conditions11,12.

Unlike time-domain and multi-sensor imaging (i.e., merging of images of the same object obtained with different imaging techniques, for example, in different wavelengths) where substantial progress was made by means of Fourier transform and wavelets1, the spatial domain yet lacks any comprehensive methodology to unify information across scales. Previous work regarding merging spatial information usually considered a single type of data, often measured by a single image, with very few material phases. Radlinski et al.13. combined SEM images and small-angle neutron scattering (SANS) data within a statistical framework to quantify the structure of porous solid in terms of continuous pore-size distribution spanning five orders of magnitude. A similar approach was used to characterize pore structure and permeability of sandstone samples using combined information from thin-sections and mercury intrusion porosimetry14. Bimodal porous structures were previously created by superimposing two truncated Gaussian random fields to simulate vuggy carbonates15. A similar methodology was also recently utilized to synthesize layered structures by connecting two separate images16. These bi-gaussian approaches can be considered a special case of the more general pluri-gaussian method popular in geostatistics. Previous works further include multiscale structures using truncated pluri-gaussian simulations17,18, however mathematically elegant, this approach utilizes only a limited number of statistical descriptors and requires tuning. Okabe and Blunt19,20 developed a general approach to a 2/3D reconstruction of the structure of porous solids via multiple-point statistics in a multi-grid framework using 2D cuts of sandstones. A slightly different Bayesian based approach was reported by Mohebi et al.21. to reconstruct double porosity features of different artificial and natural porous media samples using coarse and high resolution 2D magnetic resonance images. Continuum reconstruction of Fontainebleau sandstone was performed using pre-defined grain shapes based on multiscale digitization of sphere pack22. Other groups of researchers23,24 modelled carbonate rock structures with pore-sizes varying in size up to four orders of magnitude using SEM/thin-section and XMT data; their approach considered a primordial filter (static data from coarse scale image) decorated with crystal-filled molds of different size typical of dolomite structure. Merging two different scales of pore sizes in rocks was recently demonstrated using pore-network generation25,26,27.

The aim of this study is to develop a general framework for merging information from any number of spatial scales of any resolution into a single image, which would be simple, robust, and efficient. The methodology will be tested using structural images from shale rocks, for which data fusion will be performed using images at three spatial scales representing different spatial information on the distribution of mineral, organic matter and pore phases. Our results are expected to have practical implications in numerous disciplines, e.g., petroleum engineering, geosciences, material sciences, hydrology, soil science, biology, climatology/ecology, and remote sensing.

The multiscale data fusion method will be demonstrated on the basis of two synthetic test problems. However synthetic, these tests are aiming to solve one of the long standing challenges in petroleum geosciences, i.e. to fuse or merge structure data from macro, micro and nanoscale into a single, multiscale structure of shale rock. Shales are the source rock of shale gas and shale oil, and are known to have multiscale pore systems, composed of macroscale fractures, micropores, and nanoscale pores within gas or oil-producing organic material, referred to as kerogen. The former two scales are visible on XMT images, but are usually not connected into a percolation network28. The nanoscale pore network is organized in different interconnected clusters and extends through mudrocks where desorption and diffusion occurs within the kerogen with subsequent transport into larger-scale micropores and fractures produced by hydraulic fracturing29,30,31. Fractures may be naturally present or produced during hydraulic fracturing. XMT can efficiently resolve the 3D geometries of fractures and microscale pores within shale rocks, but cannot capture nanoporosity28. Higher-resolution two-dimensional imaging techniques such as SEM or FIB-SEM are effective to discern the pore network down to scales of several nm32. Such pore-networks find their application in the simulation of oil and gas flow during design and operation of hydrocarbon reservoirs, using so-called pore-scale models that produce effective physical properties for use in larger scale cellular Darcy-scale models28,29. To this end, the macro, micro and nanoscale pore-system have to be merged into a single network continuum with proper connections between scales.

Figure 1: Three images of shale structure representing macro, micro and nanoscale (from left to right). Below each image we report some common methods to obtain spatial information for widely differing scales. Typical values for field of view and resolution are given in parentheses.

Both examples consider merging three multiscale images with scale ratios of 16:4:1. In other words, the macroscale is four times coarser than the microscale, which in turn is four times lower in resolution than the nanoscale image. For reasons of visibility, a relatively small scale magnification factor of four was used. Although the real images used in Fig. 1 have a different scale ratio than 16:4:1, this does not preclude their use in our synthetic examples. The original image sizes were 10242 pixels for both macro and microscale and 4002 for nanoscale. Each image is first segmented into a binary image of white and black phases representing different rock phases (SI Fig. 1). We consider two cases in which unique binary phases (i.e. materials) exist at each spatial resolution. In the first case, the macroscale image represents spatial information on solids or minerals (white phase) and a mixture of microscale solids not identifiable on the macroscale and kerogen (black phase). The latter phase is resolved on the microscale image in terms of minerals (white) and poorly defined kerogen (black) phases. Finally, the nanoscale image has a sufficient resolution to discriminate the kerogen into kerogen solids (white) and kerogen porosity (black). This is a common situation where each spatial scale represents unique information on one of the material phases due to a large contrast in resolution or size of structural features. For example, kerogen nanoporosity can be resolved only on FIB-SEM/SEM images, while on XMT images kerogen is visible only as a single phase representing a mixture with other materials28.

The second case is about merging multiscale pore networks while using exactly the same input images. The pore network information includes 1) large cracks (white) and a mixture of solids and kerogen (black) on the microscale image, 2) microscale porosity (black) and kerogen phase (white) on the microscale image, and 3) similar to first case, kerogen solids (white) and kerogen porosity (black) on the nanoscale image.

Both cases represent a simplification of the real shale structure, mainly because the number of relevant phases is reduced in our study28. However, it provides us with an easily comprehensible illustration utilizing only binary information for each scale. By considering the resolution of XMT to be 1 μm and representative for the microscale, the image resolutions are fixed as follows: 4 μm for macroscale, 1 μm for microscale and 250 nm for the nanoscale image. Note that such resolutions are quite typical for routine rock imaging, as resolutions of 4–15 μm are usually obtained for standardized petrophysical cylinders used for laboratory measurements of permeability and other filtration properties10. Resolutions of 0.7–1 μm are the usual resolution for desktop scanners with cylindrical subdrilled sample with diameter of 4–10 mm28. The highest resolution is obtained with SEM images, which typically is in the range of several tens of nm32.

The backbone of the data fusion approach lies in the use of stochastic reconstruction34,35, a mathematical technique that produces a statistically similar representation of the real image and possess similar statistical properties. Previous examples have usually used the same image from one or all scales1,21,23,24, an approach that limits the size of reconstructed images or requires to use tiling (replicating similar image by connecting its copies of itself). Tiling is not useful for natural porous media as natural properties are not periodic and, thus, phases would not meet each other at the edges of tiles. Stochastic reconstructions using correlation functions is a versatile technique36 capable of dealing with any image size and with a wide range of complex structural descriptors. In addition to geometrical complexity, the image size is another important factor as it determines whether or not the image is large enough to capture all spatial features accurately.

Figure 2: Example of four times downscaled correlation functions for macroscale image. Inset shows a magnification of the beginning part of the functions.

Among many stochastic reconstruction techniques applicable to heterogeneous porous media, the more popular ones are: (i) the original method of Yeong-Torquato, combining correlation functions with simulated annealing35, (ii) the multiple-point statistics method19,20,37,38 and (iii) Gaussian random fields39. The modified Yeong-Torquato method was chosen here for reconstruction of all three scales because (i) of its ability to reproduce both isotropic and anisotropic porous media40 and (ii) the correlation functions can be rescaled to accommodate either smaller or larger scales. Spatial two-point correlation functions can take different forms (e.g. linear, cluster) and are calculated as probabilities of different image events or patterns related to a line segment of given length r41. An important feature of the correlation functions is the ability to be scaled, by which the resolution of the original correlation function can be upscaled or downscaled (Fig. 2a). Upscaling of two-point probability function dates back to the 90s34 and was routinely performed to coarsen the Guassian random field reconstructions to minimize scarce computational and memory resources of that time. Later on, both upscaling and downscaling of hypothetical two-point probability functions were used to create synthetic bimodal structures15,16. Rescaling of different types of directional correlation functions calculated from multiscale images of real porous media as part of the Yeong-Torquato stochastic reconstructions has, to the best of our knowledge, not been previously reported. Linear r interpolation is used for downscaling the macroscale image producing an increased resolution (a factor of four in Fig. 2b); a similar approach was used for upscaling the correlation function resulting in a coarsening of the resolution (SI Fig. 6). One of the main advantages of rescaling correlation functions is the ability to stochastically reconstruct images originally different in resolution into images with identical resolution, i.e., representing a super-resolution process for binary image42. This approach does not require any assumptions about unresolved spatial information (e.g., about fractal characteristics43). The latter means that smooth interpolation does not necessarily produce similar results compared to correlation functions calculated directly on the magnified image. Nevertheless, this is exactly what we want at this stage, as higher resolution details are implemented later. Downscaling and upscaling by interpolation provides reasonable accuracy for our images, as was demonstrated in a comparison between reconstruction tests on some simple images involving both rescaled correlation functions and direct image scaling (see SI for details).

The reconstruction procedure involves developing all images (three in the current example) with the same resolution: our example considers 1 μm resolution, the intermediate of three scales, to demonstrate both up and downscaling procedures. In theory, images of identical resolution can be based on any of the original scales, although the smallest-scale resolution would be preferred. Downscaling the correlation functions for the 4 μm resolution macroscale image four times produces a reconstructed image of 1 μm resolution at the microscale. To illustrate the applicability of the method to other reconstruction techniques such as multiple-point statistics, simple coarsening of the images is applied for nanoscale image before reconstruction rather than rescaling the correlation functions (the original nanoscale image resolution was coarsened four times by image rescaling, i.e., from 4002 to 1002 pixels, which produced a relatively similar reconstruction as based on a set of rescaled correlation functions with some loss of porosity, see examples and discussion in SI Fig. 7–9). Three images of equal size (40962 pixels) are generated with a resolution of 1 μm based on correlation functions for macro, micro and nanoscale images. For reasons discussed in SI (SI. 2 and SI Fig. 5), nanoscale image tiling was used, which involved a reconstructed 10242 pixel image that was multiplied four times to obtain a 40962 pixel image. By applying periodic boundary conditions during stochastic reconstruction, tiling produces smooth images without edge connection problems. Input images and resulting stochastic reconstructions are shown in Fig. 3 (separate images are also available in SI). Methodological details about stochastic reconstructions are covered in the Materials and Methods section.

Figure 3: Original segmented images (a–c) and corresponding stochastic reconstructions (d–f) used for image fusion. Original macroscale (a), microscale (b), and nanoscale (c) image. Reconstructed macroscale (d), microscale (e), and nanoscale (f) image. Insets show magnified area to compare with original image.

The overall workflow schematic for merging spatial information from macro, micro and nanoscale images is shown in SI Fig. 10. Once images of similar size and resolution have been created from rescaled correlation functions, the workflow proceeds with the actual merging step. This involves a step by step embedding of all phases according to the following sequence (based on the shale rock data for the first case). The example considers a two-step infilling of the non-mineral phase (white coloured phase in Fig. 4); in theory merging scales could also be applied to the mineral phase if subsequent higher resolution images exist with identification of different minerals. At first, we merge macro and microscale information by substituting mineral and kerogen phases (purple on Fig. 4) into the white phase on the macroscale image to obtain the combined macro-microscale image. Step two involves merging the previously obtained image with nanoscale information: the kerogen structure is substituted (orange on Fig. 4) into the white phase. Note that the white phase before each infilling step represents a phase which is not properly spatially resolved and, thus, requires information substitution from a higher-resolution imaging method. In summary, embedding is performed by overlaying three stochastic reconstruction images with decreasing scale (increasing resolution) whereby the mineral phase (blue on the macroscale and purple on the microscale) is opaque and the non-mineral phase (white) is transparent for the next higher-resolution scale (Fig. 4). For the second case of multiscale porosity, the merging procedure is similar, except that white and black phases are treated differently. The final fused image (Fig. 5) is shown in binary way to clearly highlight the multiscale pore network (shown in black) and solids (white).

Figure 4: A step by step embedding of three phases (minerals, kerogen, and pores) by fusion of three stochastic reconstructions. Phases are combined as follows: step 1 combines (i) resolved mineral (blue) and unresolved non-mineral phase (white) on macroscale image plus (ii) mixture of resolved minerals and kerogen phases (purple) and unresolved kerogen with nano-pores (white) on microscale image; step 2 combines the output of step 1 plus (iii) resolved kerogen (organic matter) phase (orange) and resolved kerogen nano-porosity (white) on nanoscale image.

Figure 5: Multiscale pore network obtained from fusing three different scales representing macrocracks, microscale and nanoscale porosity. Merging is performed in similar manner as on Fig. 4; porosity is shown in black, while white represents solids.

While the overall methodology is similar for both examples, there is a difference in the way images are merged and phases combined. The first example is specifically designed for applications where a spatial arrangement of numerous phases, materials or values is important: e.g., for 1) multiphase materials such as composites41,44, 2) soil moisture, carbon fluxes or vegetation maps11,12, 3) porosity, permeability and other hydraulic properties distributions5,7,8,45. Our second example highlights the problem of merging multiscale pore-networks that can be used for pore-scale modelling of effective properties (such as apparent and relative permeabilities10,25,26,27,28,29) or utilized for flow and transport simulations in highly heterogeneous formations3,4,46,47. An example involving both approaches can be important for modelling shale gas production rates, or for multiphase flow in mixed-wet systems where knowledge of not only the pore-network itself, but also of pore forming material is crucial to properly quantify desorption rates or menisci configurations10,30. All full size fused images from these examples are available as 40962 pixels tiff files and are included into SI.

The methodological framework for merging spatial information from different scales into a single dataset has many potential applications. The relevance of this method to earth sciences was illustrated by solving the long standing challenge of combining multiscale structural information for shale rock. The backbone of our method is based on stochastic reconstructions using rescaled correlation functions. This stochastic approach is especially useful in case high-resolution data cannot be obtained for the entire object under study, a case that exists for shales as FIB-SEM/SEM imaging is too laborious to be performed within the whole volume of the rock sample. However, for reasons of display clarity, the test problem was solved in 2D only, although the extension to 3D is straightforward48. Note that 3D problems can be solved using only a set of 2D images (cross-perpendicular 2D cuts for anisotropic 3D structures) as input data, thus highlighting another very important advantage of the stochastic approach. It is important to mention that for porous media applications such as the shale example considered here, 3D realization of our technique with real scale ratios will result in large 3D images with over a billion of voxels. With recent advancements in CFD/pore-scale modelling on the core-scale volumes49,50, our methodology is now also very appealing for upscaling of porous media properties, or for verifying other approaches such as averaging/homogenization51,52, pore-network generation25,26,27, or continuum Darcy-scale simulations28,46,47.

Other advantages of the data fusing method include: (i) provided that input data is representative for all given scales and phases, merging can be performed using limited data (i.e. a single REV size image for each embedded phase suffices), (ii) a combination of any practical number of scales and/or images can be performed, and (iii) for each of the final phases, the resulting image is statistically similar to the original images. Rescaling of correlation functions for obtaining a multiscale image has several additional benefits (based on the shale rock example): (i) it keeps phase ratios constant for any rescaling, and (ii) it does not produce “pixelized” cubical objects during downscaling (SI Fig. 7). The rescaling reconstruction can be further refined by incorporating information from higher-resolution images for cross-correlation with lower-resolution images or by performing rescaling using so-called basis functions44,53.

At this stage the abilities of the proposed method are limited mainly by 1) computational resources (as merging scales of higher resolution contrasts requires stochastic reconstructions of larger images), and 2) accuracy of the stochastic reconstructions. The first limitation can be effectively resolved by incorporating hierarchical stochastic reconstruction methods54,55, implementing some means of parallelization56 and more sophisticated pixel/voxel interchange strategies55,57, or by using tiling. Details on computing resource specifications and the computation times for the current example are provided in SI (section SI. 6). Reconstruction accuracy can be further improved by using a more comprehensive set of correlation functions to better capture image information content, as at present no universal method exists58,59. A detailed discussion on current issues related to reconstruction accuracy was recently published60 and is not reported here. It was proven15 that an intact autocorrelation (two-point probability) function for a superimposed bimodal structure is the same as for a fully resolved image. This leads us to conclude that the accuracy of our merging methodology should depend only on the accuracy of the stochastic reconstruction. The main drawback of the current realization is that cross-correlations between different phases during merging are not taken into consideration. This can be an issue if non-stationarities or transition zones are present at border regions between different phases (i.e., correlation functions for a material close to an interface is different from the bulk phase). Such information is not available in our synthetic dataset, as interfaces between embedded phases are not resolved on images (e.g., the kerogen-minerals interface is not captured on an SEM image) and requires more research. If better input data can be obtained, cross-correlations can be effectively implemented during merging procedures as a separate simulated annealing cycle. The same approach can be used for “stitching” together tiles obtained from rescaling and reconstruction of multiple input images derived from the same spatial scale.

Having numerous advantages, correlation functions are not the only stochastic reconstruction methods. As mentioned earlier, other techniques include 1) Gaussian random fields39, 2) multiple-point statistics19,20,37,38, 3) entropic descriptors61, 4) morphological algorithms23,24, and 5) process-based or grain algorithms62,63. Multiple-point methods, for instance, can be easily implemented here by scaling all input images to a fixed resolution, as was done with the nanoscale image in our test problem. Many possibilities exist to hybridize our approach. For example, for sedimentary rocks, multiple-point methods can be used for the macroscale, process-based methods can be used for the microscale and correlation functions for the nanoscale. Irrespective of the reconstruction method implemented, any number of phases can be implemented by multiphase segmentation64. For example, in our shale rock test case, microscale pores and chalcedony porosity in the mineral phase could be incorporated28. The basic idea of embedding information into different phases is shown in SI Fig. 11. This shows that however the realization of merging (or superimposing) in current examples is similar to that in the pluri-gaussian methodology, it is actually more flexible and allows for more diverse superpositioning. Tiling (or other modifications) can be also implemented at this stage to incorporate information from different images of the same phase, e.g. numerous SEM images of kerogen with different pore structure and porosity values28. Different degrees of determinism can be employed during fusing as, for example, if the coarsest information available covers the whole object under study a direct structure can be used during merging instead of stochastic reconstruction (SI Fig. 12). The number of phases and reconstruction accuracy can also be significantly increased by multiphase reconstruction methods45,55,65.

In addition to applications in the aforementioned scientific fields, the multiscale data fusion methodology can help resolve a number of other fundamental challenges, such as: 1) to study and determine REVs66 of physical parameters used (coupled with pore-scale modelling) in large-scale simulation models of flow, mass transfer and electromagnetic or seismic wave propagation, 2) to reconstruct statistically inhomogeneous random media67 by sequential reconstruction and subsequent merging of homogeneous parts, and 3) to compress voluminous images by representing them as numerous sets of correlation functions allowing reconstruction on a short notice as an efficient way of data storage and retrieval.

To stochastically reconstruct each single-scale image for a given size we utilize different sets of (rescaled) directional correlation functions40,68 computed from original binary images of real rocks. This involves the Yeong-Torquato technique35 based on simulated annealing optimization69. In this study three types of correlations functions are employed: (i) the two-point probability function S2 describing the probability that two points separated by a vector displacement r (x1, x 2) between x1 and x2 lie in the same phase, (ii) the linear function L2 describing the probability that the whole segment r lies within the given phase, and (iii) the two-point cluster function C2 describing the probability that two point separated by r lie in the same cluster. Reconstructions using S2 alone are known to be inaccurate due to the numerous degeneracy states and insufficient information content59. Simultaneous use of S2 and L2 was shown to improve reconstruction quality significantly70 and addition of C2 was effective in solving connectivity issues36,44.

It is well known that two-point probability functions do not discriminate between binary phases, this means that one cannot improve statistical information for a given structure by computing both S2(w) and S2(b) (superscript refers to the binary phase, which is either white or black). The L2 and C2 functions, unlike two-point probability functions, do discriminate between phases and can be utilized simultaneously to increase information content of the correlation function sets used for the reconstruction procedure. We calculate S2, L2 and C2 functions in two orthogonal and two diagonal directions, which were then used separately during reconstruction40.

Provided with correlation functions, we reconstruct spatial structure using the Yeong-Torquato technique, which tries to match correlation functions of a given realization with an original reference structure by pixel permutations. If a set of two-point correlation functions used in reconstruction is provided in the form of f2α(r ), where α is a type of correlation function and r is a segment of varying length, the difference between two realizations of the structure can be expressed as the sum of squared differences between sets of correlation functions:35,40,68

where f2α(r) and f2α (r) are the values of the correlation function sets for two realizations (where the former represents a reference structure while the latter represents the structure under reconstruction), wα is a weighting factor used to improve convergence. In Eq. 1 E represents the “energy” of the system, which is minimized by the simulated annealing algorithm. Initially, we create a random structure and start to change pixel positions, while checking the system’s energy according to Eq. 1. Because at the beginning of this process the characteristic sizes of phase aggregates are typically smaller than for the original image, it is reasonable to accept more permutations even if they do not reduce E. For the purpose of allowing greater flexibility at the initial stage of the inverse modelling scheme, a cooling schedule is invoked for the simulated annealing algorithm which describes the probability of accepting any permutation p in the following way:

where T is the so-called “temperature” of the system, as interpreted from the Boltzmann distribution used in the lower part of Eq. 2, and

The initial temperature T is chosen so that the probability p for ∆E ≥ 0 equals 0.548. The idea behind the cooling schedule is that simulated annealing will result in a global minimum of energy E, and would not get stuck in some local minima. From numerous trials and experience in reconstructing different test cases, we chose a slightly slow cooling schedule based on geometrical progression of the form:

where k is time step and λ is a parameter smaller but close to unity.

To accelerate simulations, we optimize recalculation of correlation functions to avoid recalculating the whole set after each pixel permutation. For this purpose we used a very efficient purpose-built optimization algorithm developed for S2 and L2 functions; a similar approach has not yet been developed for the C2 function. Although some speed-up improvements for cluster functions have been reported44, reconstructions based on small domain sizes suggests that they are still less efficient compared to S2 and L2 functions. Also, instead of random permutations, we adopted a relatively simple permutation approach following ref.70 . This involved (i) choosing a random location within a phase of interest, and (ii) choosing two random directions in which two pairs of pixels with a minimum distance in-between are selected such that they satisfy the conditions of lying in opposite phases and at the interface. Finally, further reduction in computational effort can be based on limiting the length of segments used for calculating correlation functions for the reconstruction procedure, i.e. by applying a cut-off to |r|.

Periodic boundary conditions were applied for correlation function evaluation during both the reference set evaluation and reconstruction procedure. The reconstruction procedure was terminated after 106 consecutive unsuccessful permutations. We used the annealing schedule parameter λ = 0.999999 for all reconstructions. Weighting factors wα were chosen according to the recently proposed methodology71. Rescaling was performed only for the macroscale image as explain in Fig. 2b by taking averages between adjustment points. During the reconstruction procedure each direction for each function is included separately in Eq. 1. Following sets of correlation functions were used: S2(w)−L2(w)−L2(b) for macroscale and microscale images (evaluated in four directions, resulting in a total of 12 correlation functions), and S2(w)−L2(w)− L2(b)−C2(b) for the nanoscale image (resulting in a total of 16 functions). Details about correlation function sets as well as resulting full images for each scale can be found in SI Fig. 2–5.

Macro and microscale images were binarized by the indicator kriging segmentation method72,73; the nanoscale SEM image was segmented as described in ref. 28.

1. Goshtasby, A. A. & Nikolov, S. Image fusion: Advances in the state of the art. Inf. Fusion 8(2), 114–118 (2007).

2. Prokop, M. New challenges in MDCT. Eur. Radiol. 15, E35–E45 (2005).

3. Guerriero, V. et al. Quantifying uncertainties in multi-scale studies of fractured reservoir analogues: Implemented statistical analysis of scan line data from carbonate rocks. J. Struct. Geol. 32(9), 1271–1278 (2010).

4. Strijker, G., Bertotti, G. & Luthi, S. M. Multi-scale fracture network analysis from an outcrop analogue: A case study from the Cambro-Ordovician clastic succession in Petra, Jordan. Marine Petrol. Geol. 38(1), 104–116 (2012).

5. Hamzehpour, H., Rasaei, M. R. & Sahimi, M. Development of optimal models of porous media by combining static and dynamic data: The permeability and porosity distributions. Phys. Rev. E 75(5), 056311 (2007).

6. Viceconti, M., Clapworthy, G., Testi, D., Taddei, F. & McFarlane, N. Multimodal fusion of biomedical data at different temporal and dimensional scales. Comput. Methods Programs Biomed. 102(3), 227–237 (2011).

7. Cassidy, R. et al. Combining multi-scale geophysical techniques for robust hydro-structural characterisation in catchments underlain by hard rock in post-glacial regions. J. Hydrol. 517, 715–731 (2014).

8. Jha, S. K., Comunian, A., Mariethoz, G. & Kelly, B. F. J. Parameterization of training images for aquifer 3-D facies modelling integrating geological interpretations and statistical inference. Water Resour. Res. 50(10), 7731–7749 (2014).

9. Khalili, A. D. et al. Permeability upscaling for carbonates from the pore scale by use of multiscale X-Ray-CT images. SPE Reservoir Eval. Eng. 16(4), 353–368 (2013).

10. Blunt, M. J. et al. Pore-scale imaging and modelling. Adv. Water Resour. 51, 197–216 (2013).

11. Wang, Y. P., Trudinger, C. M. & Enting, I. G. A review of applications of model-data fusion to studies of terrestrial carbon fluxes at different scales. Agric. For. Meteorol. 149(11), 1829–1842 (2009).

12. Bradley, E. S., Toomey, M. P., Still, C. J. & Roberts, D. A. Multi-scale sensor fusion with an online application: integrating GOES, MODIS, and webcam imagery for environmental monitoring. IEEE J. Sel. Top. Appl. Earth Observ. Remote Sens. 3(4), 497–506 (2010).

13. Radlinski, A. P. et al. Angstrom-to-millimeter characterization of sedimentary rock microstructure. J. Colloid Interface Sci. 274, 607–612 (2004).

14. Amirtharaj, E. S., Ioannidis, M. A., Parker, B. & Tsakiroglou, C. D. Statistical synthesis of imaging and porosimetry data for the characterization of microstructure and transport properties of sandstones. Transp. Porous Media 86(1), 135–154 (2011).

15. Moctezuma-Berthier, A., Vizika, O. & Adler, P. M. Macroscopic conductivity of vugular porous media. Transp. Porous Media 49(3), 313–332 (2002).

16. Hyman, J. D. & Winter, C. L. Stochastic generation of explicit pore structures by thresholding Gaussian random fields. J. Comp. Phys. 227, 16–31 (2014).

17. Dowd, P. A., Pardo-Iguzquiza, E. & Xu, C. Plurigau: a computer program for simulating spatial facies using the truncated plurigaussian method. Comput. Geosci. 29(2), 123–141 (2003).

18. Mariethoz, G., Renard, P., Cornaton, F. & Jaquet, O. Truncated plurigaussian simulations to characterize aquifer heterogeneity. Ground Water 47(1), 13–24 (2009).

19. Okabe, H. & Blunt, M. J. Prediction of permeability for porous media reconstructed using multiple-point statistics. Phys. Rev. E 70(6), 066135 (2004).

20. Okabe, H. & Blunt, M. J. Pore space reconstruction using multiple-point statistics. J. Petrol. Sci. Eng. 26(1–2), 121–137 (2005).

21. Mohebi, A., Fieguth, P. & Ioannidis, M. A. Statistical fusion of two-scale images of porous media. Adv. Water Resour. 32(11), 1567–1579 (2009).

22. Latief, F. D. E., Biswal, B., Fauzi, U. & Hilfer, R. Continuum reconstruction of the pore scale microstructure for Fontainebleau sandstone. Physica A 389(8), 1607–1618 (2010).

23. Biswal, B., Oren, P. E., Held, R. J., Bakke, S. & Hilfer, R. Stochastic multiscale model for carbonate rocks. Phys. Rev. E 75(6), 061303 (2007).

24. Roth, S. et al. Continuum-based rock model of a reservoir dolostone with four orders of magnitude in pore sizes. Am. Assoc. Pet. Geol. Bull. 95(6), 925–940 (2011).

25. Jiang, Z, van Dijke, M.I.J., Sorbie, K.S. & Couples, G.D. Representation of multiscale heterogeneity via multiscale pore networks. Water Resour. Res. 49, 5437–5449 (2013).

26. Bauer, D. et al. From computed microtomography images to resistivity index calculations of heterogeneous carbonates using a dual-porosity pore-network approach: Influence of percolation on the electrical transport properties. Phys. Rev. E 84, 011133 (2011).

27. Bultreys, T., Van Hoorebeke, L. & Cnudde, V. Multi-scale, micro-computed tomography-based pore network models to simulate drainage in heterogeneous rocks. Adv. Water Resour. 78, 36–49.

28. Gerke, K. M. et al. Determining physical properties of unconventional reservoir rocks: from laboratory methods to pore-scale modeling. SPE 167058 Technical paper, doi: 10.2118/167058-MS (2013).

29. Mehmani, A., Prodanović, M. & Javadpour, F. Multiscale, multiphysics network modeling of shale matrix gas flows. Transp. Porous Media 99(2), 377–390 (2013).

30. Mehmani, A. & Prodanović, M. The application of sorption hysteresis in nano-petrophysics using multiscale multiphysics network models. Int. J. Coal Geol. 128, 96–108 (2014).

31. Patzek, T. W., Male, F. & Marder, M. Gas production in the Barnett Shale obeys a simple scaling theory. Proc. Natl. Acad. Sci. USA 110, 19731–19736 (2013).

32. Loucks, R. G., Reed, R. M., Ruppel, S. C. & Jarvie, D. M. Morphology, genesis, and distribution of nanometer-scale pores in siliceous mudstones of the Mississippian Barnett shale. J. Sediment Res. 79(11–12), 848–861 (2009).

33. lmishek, G. F. Petroleum geology and resources of the West Siberian Basin, Russia”. US Geological Survey Bulletin 2201-G, 49 (2003).

34. Adler, P. M., Jaquin, C. G. & Quiblier, J. A. Flow in simulated porous media. Int. J. Multiphase Flow 16(4), 691–712 (1990).

35. Yeong, C. L. Y. & Torquato, S. Reconstructing random media. Phys. Rev. E 57, 495–506 (1998).

36. Jiao, Y., Stillinger, F. H. & Torquato, S. A superior descriptor of random textures and its predictive capacity. Proc. Natl. Acad. Sci. USA 106, 17634 (2009).

37. Mariethoz, G., Renard, P. & Straubhaar, J. The direct sampling method to perform multiple-point geostatistical simulations. Water Resour. Res. 46, W11536 (2010).

38. Tahmasebi, P. & Sahimi, M. Cross-correlation function for accurate reconstruction of heterogeneous media. Phys. Rev. Lett. 110(7), 078002 (2013).

39. Hesse, F., Prykhodko, V., Schluter, S. & Attinger, S. Generating random fields with a truncated power-law variogram: A comparison of several numerical methods. Environ. Model. Softw. 55, 32–48 (2014).

40. Gerke, K. M., Karsanina, M. V., Vasilyev, R. V. & Mallants, D. Improving pattern reconstruction using correlation functions computed in directions. Europhys. Lett. 106(6), 66002 (2014).

41. Torquato, S. Random heterogeneous materials: microstructure and macroscopic properties (Springer-Verlag, New York, 2002).

42. Nasrollahi, K. & Moeslund, T. B. Super-resolution: a comprehensive survey. Mach. Vis. Appl. 25(6), 1423–1468 (2014).

43. Mariethoz, G., Renard, P. & Straubhaar, J. Extrapolating the fractal characteristics of an image using scale-invariant multiple- point statistics. Math. Geosci. 43(7), 783–797 (2011).

44. Jiao, Y. & Chawla, N. Modeling and characterizing anisotropic inclusion orientation in heterogeneous material via directional cluster functions and stochastic microstructure reconstruction. J. Appl. Phys. 115(9), 093511 (2014).

45. Schluter, S. & Vogel, H. J. On the reconstruction of structural and functional properties in random heterogeneous media. Adv. Water Resour. 34(2), 314–325 (2011).

46. Matthäi, S. K., Nick, H. M., Pain, C. & Neuweiler, I. Simulation of solute transport through fractured rock: a higher-order accurate finite-element finite-volume method permitting large time steps. Transp. Porous Media 83(2), 289–318 (2010).

47. Wolff, M., Flemisch, B. & Helmig, R. An adaptive multiscale approach for modeling two-phase flow in porous media including capillary pressure. Water Resour. Res. 49(12), 8139–8159 (2013).

48. Yeong, C. L. Y. & Torquato, S. Reconstructing random media: II. Three-dimensional media from two-dimensional cuts. Phys. Rev. E 58(1), 224–233 (1998).

49. Khirevich, S., Ginzburg, I. & Tallarek, U. Coarse- and fine-grid numerical behavior of MRT/TRT lattice-Boltzmann schemes in regular and random sphere packings. J. Comp. Phys. 281, 708–742 (2015).

50. Aghaei, A. & Piri, M. Direct pore-to-core up-scaling of displacement processes: Dynamic pore network modeling and experimentation. J. Hydrol. 522, 488–509 (2015).

51. Garibotti, C. R. & Peszynska, M. Upscaling non-Darcy flow. Transp. Porous Media 80(3), 401–430 (2009).

52. Wargo, E. A., Hanna, A. C., Cecen, A., Kalidindi, S. R. & Kumbur, E. C. Selection of representative volume elements for pore- scale analysis of transport in fuel cell materials. J. Power Sources 197, 168–179 (2012).

53. Jiao, Y., Stillinger, F. H. & Torquato, S. Modeling heterogeneous materials via two-point correlation functions: Basic principles. Phys. Rev. E 76, 031110 (2007).

54. Campaigne, W. R. & Fieguth, P. W. Frozen-State Hierarchical Annealing. IEEE Trans. Image Process. 22(4), 1484–1495 (2013).

55. Chen, D. D., Teng, Q. Z., He, X. H., Xu, Z. & Li, Z. J. Stable-phase method for hierarchical annealing in the reconstruction of porous media images. Phys. Rev. E 89(1), 013305 (2014).

56. Chu, K. W., Deng, Y. F. & Reinitz, J. Parallel simulated annealing by mixing of states. J. Comput. Phys. 148(2), 646–662 (1999).

57. Pant, L. M., Mitra, S. K. & Secanell, M. Stochastic reconstruction using multiple correlation functions with different-phase- neighbor-based pixel selection. Phys. Rev. E 90, 023306 (2014).

58. Gommes, C., Jiao, Y. & Torquato, S. Density of states for a specified correlation function and the energy landscape. Phys. Rev. Lett. 108(8), 080601 (2012).

59. Gommes, C. J., Jiao, Y. & Torquato, S. Microstructural degeneracy associated with a two-point correlation function and its information content. Phys. Rev. E 85(5), 051140 (2012).

60. Karsanina, M. V., Gerke, K. M., Skvortsova, E. B. & Mallants, D. Universal spatial correlation functions for describing and reconstructing soil microstructure. PLoS ONE 10(5), e0126515 (2015).

61. Piasecki, R. Microstructure reconstruction using entropic descriptors. Proc. R. Soc. Lond. A Math. Phys. Sci. 467(2127), 806–820 (2011).

62. Øren, P. & Bakke, S. Reconstruction of Berea sandstone and pore-scale modeling of wettability effects. J. Petrol. Sci. Eng. 39, 177–199 (2002).

63. Thovert, J. F. & Adler, P. M. Grain reconstruction of porous media: Application to a Bentheim sandstone. Phys Rev E 83(5), 056116 (2011).

64. Schluter, S., Sheppard, A., Brown, K. & Wildenschild, D. Image processing of multiphase images obtained via X- ray microtomography: A review. Water Resour. Res. 50(4), 3615–3639 (2014).

65. Rozman, M. G. & Utz, M. Efficient reconstruction of multiphase morphologies from correlation functions. Phys. Rev. E 63(6), 066701 (2001).

66. Zeman, J. & Sejnoha, M. From random microstructures to representative volume elements. Modelling and Simulation in Materials Science and Engineering 15(4), S325–S335 (2007).

67. Tahmasebi, P. & Sahimi, M. Reconstruction of nonstationary disordered materials and media: Watershed transform and cross-correlation function. Phys. Rev. E 91(3), 032401 (2015).

68. Kumar, H., Briant, C. L. & Curtin, W. A. Using microstructure reconstruction to model mechanical behavior in complex microstructures. Mech. Mater. 38, 818–832 (2006).

69. Kirkpatrick, S., Gelatt, C. D. & Vecchi, M. P. Optimization by simulated annealing. Science 220(4598), 671–680 (1983).

70. Čapek, P., Hejtmánek, V., Kolafa, J. & Brabec, L. Transport properties of stochastically reconstructed porous media with improved pore connectivity. Transp. Porous Media 88(1), 87–106 (2011).

71. Gerke, K. M. & Karsanina, M. V. Improving stochastic reconstructions by weighting correlation functions in an objective function. Europhys. Lett. 111, 56002 (2015).

72. Oh, W. & Lindquist, B. Image thresholding by indicator kriging. IEEE Trans. Pattern. Anal. Mach. Intell. 21, 590–602 (1999).

73. Houston, A. N., Otten, W., Baveye, P. C. & Hapca, S. Adaptive-window indicator kriging: A thresholding method for computed tomography images of porous media. Comput. Geosci. 54, 239–248 (2013).

We are grateful to Brendan Speet (CSIRO Land and Water) for his help with computing resources. We also thank Prof. Vladimir Voevodin for providing computational resources of Moscow State University Computational Center. This work was partially supported by RSF grant 14-17-00658 (stochastic reconstructions using directional correlation functions), RFBR grant 15-34-20989 мол_а_вед (data fusion), and RFBR grants 13-04-00409-a and 13-05-01176-a. We are grateful to Grégoire Mariéthoz for fruitful discussions. Special thanks to Dina Gilyazetdinova, Dmitry Korost, Natalya Balushkina and Georgy Kalmykov for all shale rock images, and to Timofey Sizonenko for segmentation help. Comments on an earlier version of this manuscript have been obtained from Tobias Mueller (CSIRO Energy) and Ben Clennell (CSIRO Energy).

K.M.G. conceived the idea. K.M.G. and M.V.K. performed simulations. K.M.G., M.V.K. and D.M. wrote the manuscript.

Competing financial interests: The authors declare no competing financial interests.

How to cite this article: Gerke, K. M. et al. Universal Stochastic Multiscale Image Fusion: An Example

Application for Shale Rock. Sci. Rep. 5, 15880; doi: 10.1038/srep15880 (2015).

1CSIRO Land and Water, Glen Osmond, PB2, SA 5064, Australia. 2The University of Melbourne, Department of Infrastructure Engineering, Parkville, VIC, 3010, Australia. 3Institute of Geosphere Dynamics of the Russian Academy of Sciences, Leninsky prosp. 38/1, Moscow, 119334, Russia. 4Institute of Physics of the Earth of Russian Academy of Sciences, Bolshaya Gruzinskaya 10, Moscow, 107031, Russia. *These authors contributed equally to this work. Correspondence and requests for materials should be addressed to K.G. (email: kirill.gerke@unimelb. edu.au)

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

BOOKMARKS

AllAboutShale.com