Data-constrained characterization of sandstone microstructures with multi-energy X-ray CT

A data-constrained non-linear optimization approach has been developed to characterize microscopic distributions of mineral phases and pores in a sandstone sample using X-ray CT data sets acquired at 35keV and 45keV beam energies as constraints. The approach minimizes discrepancy between the expected and measured linear absorption coefficients and maximizes Boltzmann distribution probability. It enables integration of both the 3D X-ray CT data-constraints and global level information, and leads to more accurate predictions of microscopic 3D compositional distributions in material samples. Permeability simulations and comparisons with experimentally measured porosity indicate that DCM characterisation agrees reasonably with experimental observations. However, segmentation of CT images leads to under-estimation of porosity and permeability.


Introduction
In oil and gas reservoir rock characterization, and in materials development and processing, it is often the case that one needs to optimize a given set of material performance properties.Material performance properties are closely related to their microstructures.Considerable experimental, theoretical and numerical efforts have been devoted to 3D microscopic characterization [1][2][3][4][5][6][7].X-ray CT has been widely used as a sample non-destructive method for microstructures characterization.However, segmentation of X-ray CT images is not always adequate to discriminate material compositions, as different combinations of material phases may exhibit similar X-ray absorption properties.
A level of success has been achieved in microstructures characterization using a data-constrained modeling (DCM) linear optimization approach with multi-energy (spectrum) X-ray CT [8][9][10][11][12][13][14][15][16][17].However, it is relatively sensitive to experimental noise, and in general not suitable for characterization of materials consisting of more than three compositions.In this article, both multienergy X-ray CT data and microscopic correlations among material compositions are incorporated into the DCM model.The correlations are formulated as phenomenological self energies for each composition as well as intra-and inter-voxel interaction energies between different compositions.The overall formulation of the problem resembles one of statistical mechanics [18].It is capable of treating partial occupancy of voxels and can in principle incorporate a wide range of types of data and sources of a priori information about the sample.The methodology is generic and is applicable to a wide range of science and engineering disciplines.

DCM nonlinear optimization formulation
The model is defined on a simple-cubic grid with a cubic voxel (volumetric pixel) on each grid point.
For a voxel at n ( ) where N is the total number of voxels in the system, the model attempts to minimize the following objective function: where L is the number of X-ray CT data sets acquired from the same sample using L different beam energies.This is equivalent to minimizing the discrepancy between the expected and the measured linear absorption coefficients and maximizing Boltzmann distribution probability [18][19] For each X-ray spectrum (beam energy) l , the relative CT error is expressed as . The optimization is achieved by adjusting the volume fraction variables ) for each material composition m , where M is the total number of compositions in the system, subject to the constraints: The dimensionless phenomenological interaction energy is expressed as  Figure 1 shows that calcite exists both as coatings on the quartz grains and as clusters between quartz grains.This is consistent with SEM observations [19].In DCM generated images, calcite pixels have varying intensities.This indicates the coexistence of calcite and void.In other words, calcite exists in varying degrees of porosity in the sample which is consistent with experimental observations.This partial occupancy of voxels cannot be obtained by image segmentation.

Permeability simulations
A random walk method is used to predict permeability defined by Darcy's law [20].In general, permeability is a function of compositional volume fractions and internal fine structure.For the predefined sandstone sample, compositional volume fractions have been derived from the DCM nonlinear model or least-square segmentation model.However, the fine structure below the voxel size is unknown.An analytic formula derived in [20] is used to calculate voxel based permeability under an isotropic structure constraint, i.e. compositions mixed randomly.Our numerical study in this section will use an isotropic structure for the effective permeability simulations.Figure 2 shows the permeability simulation results based on different pore connectivity constraints.For each constraint, simulated permeability is plotted against the inverse of simulation physical scale.Circles on the vertical axis show converging values of permeabilities for different pore configurations.Since quartz is not permeable and the permeability of calcite is insignificant compared to void, only void permeability is required to calculate voxel-based permeability.In our numerical simulations, the void phase has been assigned a dimensionless permeability value of 100, and quartz and calcite have been assigned a value of 0. The permeability value of a voxel is calculated based on its compositional volume fractions [20].The top and bottom plots in the figure represent the simulations based on the theoretical maximum and minimum to define an upper and a lower bound for the effective permeabilities.The middle three plots correspond to three pore structure configurations: isotropic mixture of compositions, from DCM non-linear optimization and from least-square segmentation.Conclusions to draw from Figure 2 are firstly, the simulated permeabilities based on different pore structure constraints are bounded by the theoretical range, thereby verifying the random walk algorithm.Secondly, the least-square segmentation method leads to a reduction in pore connectivity, which results in a significant underestimation of effective permeability.Thirdly, DCM non-linear optimization improves pore network activity.The simulated permeability from the DCM output shows a quantitative agreement with the permeability simulated by applying an isotropic structure constraint.As can be seen, both limits are very close to each other.

Summary
Phenomenological intra-and inter-voxel interactions among material compositions as well as self-energies are incorporated into the data-constrained model (DCM) for characterization of microstructures of various materials.It is possible to incorporate prior knowledge and other measurements.The model is formulated as a constrained non-linear optimization problem.Numerical results indicate that the DCM non-linear optimization approach produces a more robust and accurate prediction of compositional microstructures, and leads to more accurate modeling of permeability properties.
The X-ray CT images and DCM reconstructed compositional maps are available from the CSIRO Data Access Portal [21].This work is sponsored by CSIRO Computational & Simulation Science Transformational Capability Platform and Shanxi Province One-Hundred Person project.
. The parameter ζ ( 0 ≥ ) is introduced to approximate relative error weight distribution for the CT reconstructed absorption coefficients.Its default value is 1 = ζ .
, ( l m µ is the total linear absorption coefficient for material m , ) ( ˆl n µ is the CT reconstructed linear absorption coefficient, and

S
denotes the self-energy (chemical potential) for composition m , energy between neighboring voxel compositions 1 m and 2 m at distance k , where 0 = k ( 0 > k ) denotes intra-(inter-) voxel.The maximum interaction range is denoted by K .The number of neighboring voxels with neighboring range k is denoted by ) (k N and the j th neighboring vector with a neighboring range k is denoted by ) (k j n .The parameter η with a default value of 5 .0 = η is used to approximate the volume to interface relation.Multi-spectrum X-ray CT with a synthetic CIPS (Calcite In-situ Precipitation System) sandstone sample has been carried out at beam energies 35 and 45keV[15].DCM non-linear optimization has been applied to the data sets which represent a region of typical CT image, DCM non-linear reconstructed and least-square segmented compositional maps are shown in Figure1.DCM non-linear reconstruction shows an average porosity of 19.6% which compares favorably with the He gas measured value of 19.5%.The model uses parameters: -energy CT least-square segmentation, which is similar to threshold segmentation for single energy CT, gives a much lower overall porosity of 11.3%.To minimize error near sample boundary and image artifacts slices 1-50 and 215-250 have been excluded in porosity estimations.

Figure 1 .
Figure 1.X-ray CT image and compositional maps at slice 140.(a) X-ray CT image at monochromatic beam energy 35keV.(b) X-ray CT image at monochromatic beam energy 45keV.(c) Microscopic compositional distributions from DCM non-linear optimization.(d) Microscopic compositional distributions from least-square segmentation.In (c) and (d), quartz, calcite and void are represented as blue, red and white respectively.Volume fractions are represented by their colour intensities.

Figure 2 .
Figure 2. Comparison of the random walk simulation results from different pore structure constraints with given  (void) = 100.Simulations based on the theoretical maximum and minimum are plotted on the top and bottom respectively.The converging point for each case is represented by a circle on the vertical axis.The simulation result from the DCM non-linear output correlates well with the isotropic microstructure constraint.