2305.06519.pdf

Homogenizing elastic properties of large digital rock images by combining CNN

with hierarchical homogenization method

Rasool Ahmad

∗

, Mingliang Liu

, Michael Ortiz

, Tapan Mukerji

, Wei Cai

Micro and Nano Mechanics Group, Department of Mechanical Engineering, Stanford University, CA 94305, USA

Department of Energy Science and Engineering, Stanford University, CA 94305, USA

Division of Engineering and Applied Science, California Institute of Technology, CA 91125 Pasadena, USA

Abstract

Determining e

ective elastic properties of rocks from their pore-scale digital images is a key goal of digital rock

physics (DRP). Direct numerical simulation (DNS) of elastic behavior, however, incurs high computational cost;

and surrogate machine learning (ML) model, particularly convolutional neural network (CNN), show promises to

accelerate homogenization process. 3D CNN models, however, are unable to handle large images due to memory

issues. To address this challenge, we propose a novel method that combines 3D CNN with hierarchical homogenization

method (HHM). The surrogate 3D CNN model homogenizes only small subimages, and a DNS is used to homogenize

the intermediate image obtained by assembling small subimages. The 3D CNN model is designed to output the

homogenized elastic constants within the Hashin-Shtrikman (HS) bounds of the input images. The 3D CNN model

is first trained on data comprising equal proportions of five sandstone (quartz mineralogy) images, and, subsequently,

fine-tuned for specific rocks using transfer learning. The proposed method is applied to homogenize the rock images of

size 300

300

300 and 600

600

600 voxels, and the predicted homogenized elastic moduli are shown to agree with

that obtained from the brute-force DNS. The transferability of the trained 3D CNN model (using transfer learning) is

further demonstrated by predicting the homogenized elastic moduli of a limestone rock with calcite mineralogy. The

surrogate 3D CNN model in combination with the HHM is thus shown to be a promising tool for the homogenization

of large 3D digital rock images and other random media.

Keywords:

Digital rock physics, Deep learning, Elastic properties, Hierarchical homogenization method, 3D

Convolutional neural networks

1. Introduction

Digital Rock Physics (DRP) aims to complement

replace the expensive laboratory experiments to characterize the

relevant properties of reservoirs directly from the pore-scale digital images of the constituent rocks [1, 2, 3, 4, 5, 6].

The DRP workflow starts from the acquisition of the three-dimensional digital images of rocks using modern imaging

techniques such as micro computed-tomography (micro-CT). The digital image is then segmented into its di

erent

constituents and pores using an image processing tool. A physical simulation is then performed numerically on the

segmented digital images to obtain the desired e

ective properties of the rock. The e

ective elastic properties, e.g.

bulk and shear moduli, are important properties of rocks to fully characterize the reservoirs and is the focus of the

present work.

Rocks are quintessential examples of random media in which microstructures (distribution of constituents materials

and pores) span multiple length scale. Resolution of the finer length-scale features of the microstructure is necessary

to obtain accurate results from the numerical simulations. On the other hand, a large field-of-view (FOV) is needed

to sample statistically representative volume to obtain well-converged results. Satisfying both conditions is becoming

increasingly possibly with the advancement in imaging technology. Furthermore, machine learning-based fusion of

high resolution and small FOV scanning electron microscopy (SEM) images with low resolution and large FOV micro-

CT images yields huge images (with multibillion voxels) with high resolution and large FOV [7]. However, performing

∗

corresponding author

Email addresses:

rasool@stanford.edu

( Rasool Ahmad ),

mliu9@stanford.edu

( Mingliang Liu ),

ortiz@aero.caltech.edu

(Michael Ortiz),

mukerji@stanford.edu

(Tapan Mukerji),

caiwei@stanford.edu

(Wei Cai)

Preprint submitted to

May 12, 2023

arXiv:2305.06519v1 [physics.geo-ph] 11 May 2023

brute-force numerical simulations on such huge images to determine e

ective properties remains out of reach due to

requirement of prohibitively expensive requirement of computational resources

−

in terms of both computation time

and memory to store images and intermediate computation steps.

Domain decomposition-based methods are commonly utilized to e

ffi

ciently carry out the numerical simulations on

huge computational domains (digital images in the present context). The basic idea of domain decomposition methods

is to first divide the huge image under investigation into multiple small computationally tractable subimages. The

numerical simulations are performed on individual subimages to obtain their e

ective properties. The final step is to

integrate the homogenized properties of the individual subimages to determine the e

ective property of the original

huge image. One way to realize this last integration is by simply computing the simple arithmetic

geometric

harmonic

mean of the properties of the individual subimages. Hierarchical homogenization method (HHM) is another princi-

pled renormalization method-inspired procedure to perform the last step of the integration to find accurate e

ective

properties of the huge domain. In HHM, the last integration step is carried out by assembling the subimages with

corresponding e

ective properties and performing the numerical simulation on much smaller partially homogenized

domain. Recently, Ahmad et al. [8] systematically analyzed the errors incurred by the HHM, and provide a procedure

to obtain accurate and unbiased e

ective elastic moduli of digital rocks.

In domain decomposition methods such as HHM, the computational cost associated with the last integration step

is negligible. The majority of the computational cost in HHM is incurred during the homogenization of the subimages

by performing relatively expensive numerical simulations. Thus, replacing the expensive numerical simulations with

cheaper deep learning-based surrogate models to predict the e

ective elastic moduli of subimages would significantly

increase the overall computational e

ffi

ciency of the HHM, and is the main focus of the present work.

Deep learning-based surrogate models have found wide-ranging applications in predicting the physical properties

of digital rocks. The strength of deep learning methods originates from their ability to learn important features from

raw data to predict the target properties. Convolutional Neural Networks (CNN) is at the core of most prominent

surrogate deep leaning models for DRP applications which operate directly on the digital images to predict the desired

properties. For instance, CNN-based models have been used for image segmentation [9, 10, 11, 12], predicting the flow

properties and permeability [13, 14, 15, 16, 17, 18, 19], predicting mechanical properties [11, 20, 21, 22, 23, 24, 25],

enhancing the resolution of digital images [7, 26, 27, 28], image reconstruction [29, 30], classification [31].

The application of CNN to large three-dimensional digital images is severely limited by the excessive memory

requirements to fit the images on graphical processing units (GPUs). Large images restrict the batch size and lead

to ine

ffi

cient training of the model. To mitigate this memory issue, Kashefi and Mukerji [32] use PointNet instead of

CNN to predict the permeability from rock images. PointNet operates on the point cloud representation of the interface

between pore and mineral phases and thus requires much lower memory than the full image. Furthermore, Santos et al.

[17] propose a multiscale CNN to predict permeability of rock images with more than 512

voxels. However, the large

image size results in the smaller number of training data which can lead to overfitting and poor generalizability of the

trained model.

In this work, we propose a novel hybrid approach, CNN-HHM, which combines deep learning model and direct

numerical simulation (DNS) in the framework of HHM to e

ffi

ciently determine the e

ective elastic moduli of large

digital rock images. The framework starts with the partitioning of a large rock image of size

voxels into

(

)

disjoint smaller subimages of size

voxels. A surrogate 3D-CNN model is trained to predict the e

ective

isotropic elastic moduli (bulk and shear) of subimages which are numerous to train a deep learning model. To build

the training data, the e

ective elastic moduli of the subimages are determined by solving the elasticity problem using

an e

ffi

cient FFT-based numerical solver. The main aspect of our 3D-CNN implementation is that it does not directly

predict the bulk and shear moduli from the segmented subimages; instead the 3D-CNN model outputs two scalars

and

having values between 0 and 1. The two numbers

and

specify the e

ective bulk and shear moduli with

respect to the respective lower and upper Hashin-Strikman bounds. This particular feature of the model results in high

accuracy and excellent generalizability of the trained model. The surrogate CNN model is pretrained on a mixed data

set containing equal proportion of five sandstone rock (quartz mineralogy) images of size 75

75 voxels. We,

subsequently, use transfer learning to fine tune the surrogate CNN model for each rock separately.

The trained 3D-CNN model is then integrated into the HHM scheme to determine the e

ective elastic moduli

of large rocks. In particular, the e

ective moduli of the subimages are obtained using the trained 3D-CNN model

instead of DNS. The substitution of computationally costly numerical simulation with the cheaper surrogate 3D-CNN

model leads to significant increase in the computational e

ffi

ciency of the overall HHM scheme. The second and

the final homogenization of much smaller intermediate assembled image is performed by direct physical simulation

which has negligible computational cost. The proposed framework is shown to predict the e

ective elastic moduli of

digital rock images of size 300

300

300 and 600

600

600 voxels in excellent agreement with DNS values.

The CNN-HHM approach with transfer learning is then shown to predict the elastic moduli of limestone (calcite

mineralogy) with reasonable accuracy. Thus, the proposed CNN-HHM approach is shown to be a promising tool for

the e

ffi

cient determination of the e

ective elastic moduli of large digital rock images with varies microstructure and

mineral composition.

The rest of the manuscript is organized as follows. In Section 2, we present the details of the hierarchical homog-

enization method (HHM) and the surrogate 3D-CNN model which forms the component of the CNN-HHM approach.

In Section 3, we present the results of the proposed CNN-HHM approach. Section 4summarizes and discusses the

main findings of the present work.

2. Methods

In this section, we first briefly describe the hierarchical homogenization method (HHM) closely following Ahmad

et al. [8]. We then present the details of the surrogate 3D-CNN model including the architecture, the loss function used

to train the model, and the generation of the training, validation and testing data.

2.1. Hierarchical Homogenization Method (HHM)

Hierarchical homogenization method (HHM) is a renormalization inspired approximate method to determine the

homogenized properties of large computational domains. This method is especially helpful for digital rock physics

application, as requirement of high resolution and large field of view result in large images that cannot be handled

by brute force direct numerical simulations. Instead of solving the elasticity PDEs in the large image, HHM adopts

the renormalization inspired approach where the large computational domain is scaled down by successive coarse-

graining.

homogenization

of red subimages

homogenization

of the assembed

image using DNS

N/n

Original 'big' image partitioned

into 'smaller' red subimages

Assembled image after one

step of homogenization

Fully homogenized image

assembling into

intermediate image

Homogenization of (N/n)

subimages of size n

DNS / CNN

Figure 1: Schematic of the HHM scheme. The large image is partitioned into smaller red subimages. The red subimages are homogenized separately

by direct numerical simulation (DNS) or a surrogate convolutional network (CNN). The homogenized subimages are assembled to obtain a partially

homogenized image. The assembled partially homogenized image is again homogenized by directly solving the elasticity PDEs to find the final

homogenized elastic constants.

HHM starts with the partitioning of a 3D large image of size (number of voxels)

into (

)

smaller

subimages of size

as schematically depicted in red in the first panel of Figure 1. Each voxel of the image

represents one of the two phases (pore and mineral in our case). The two phases are considered to be homogeneous

isotropic elastic materials with their own elastic constants. We then homogenize the red subimages and replace them

with one voxel of corresponding homogenized elastic constants. The subimages are homogenized by solving the

elasticity PDEs using periodic boundary conditions. In the next step, we assemble the homogenized subimges to

obtain a partially homogenized version of the original large image, as shown in the third panel of Figure 1. The partially

homogenized image is much smaller ((

)

voxels) in size compared to original large image, and is amenable to brute

direct numerical simulations. Furthermore, the partially homogenized image contains the low-frequency information

of the original image and the high-frequency local information are summarized into the homogenized properties of

subimages. The assembled image is again homogenized by solving the elasticity PDEs using periodic boundary

conditions to find the final homogenized isotropic elastic constants.

In the original HHM scheme, the individual subimages are extracted from the large image and homogenized by

direct numerical simulation (DNS) as depicted in Figure 2. The subimages are subjected to a periodic boundary

condition and the equilibrium stress and strain field is determined by solving the elasticity PDEs using FFT-based

elasticity solver. The e

ective elastic sti

ness matrix is then determined by relating the average stress and strain in

the subimages. To find all the components of the e

ective sti

ness matrix, the elasticity PDEs are solved for six

independent boundary conditions. For more details, reader are directed to Ahmad et al. [8].

periodic boundary condition

homogenization

DNS/CNN

Figure 2: Schematic of the homogenization of a single subimage. The subimage is composed of two phases (pore and mineral). The segmented

subimage is homogenized either by solving the elasticity PDEs (DNS) under periodic boundary conditions or using a surrogate CNN model.

Thus, HHM enables the homogenization of a large image without solving the elasticity equations in the large

image itself. Since most of the time of the HHM is spent to homogenize (

)

subimages, a fast surrogate CNN

model to obtain the homogenized elastic constant of smaller subimages will significantly accelerate the whole HHM

scheme.

2.2. Details of the surrogate three-dimensional CNN model

As depicted in Figure 2, the smaller red subimages are extracted from the large image and processed by a surrogate

3D CNN model. In this work, our 3D CNN model is designed to take images of size 75

75 as input and output

the e

ective bulk and shear moduli of the subimage. The input images are already segmented into two phases of pore

and mineral.

The architecture of the surrogate CNN model is depicted schematically in Figure 3; and the details of the various

layers involved in the model are presented in Table 1. The input image is first passed through a series of 3D convolution

layers with batch normalization, average pooling and PRelu non-linear activation layers. The output of the last 3D

convolution layer is then flattened and passed through a series of fully connected linear layers with batch normalization

and PRelu non-linear activation layers. The output of the final linear layer is passed through 0

5(tanh(

)

1) function

to constrain the values between 0 and 1. Thus, the final output of the surrogate CNN model is an array of two numbers

∈

1) and

∈

1) which indicate the e

ective moduli location between their respective Hashin-Shtrikman

bounds. Specifically, the e

ective bulk and shear moduli are obtained by using the following equations

lower

−

)

upper

lower

−

)

upper

(1)

where

lower

and

upper

are the lower and upper Hashin-Shtikman bounds of the bulk modulus, and

lower

and

upper

are the lower and upper Hashin-Shtikman bounds of the shear modulus. The use of the Hashin-Shtrikman bounds to

inform the surrogate CNN model improves the accuracy of the prediction of the e

ective elastic moduli, and enables

the transferability of the model to di

erent rock types as will be demonstrated later in Section []. Furthermore, the

Hashin-Shtrikman bounds are cheap themselves to compute from the information of the volume fraction of the two

phases (porosity in the case of rock) in the subimages and their respective bulk and shear moduli [33], as

upper

mineral

−

(

pore

−

mineral

)

−

(

mineral

)

lower

pore

(

mineral

−

pore

)

−

(

−

)

(

pore

)

upper

mineral

−

(

pore

−

mineral

)

−

(

mineral

)

mineral

(

mineral

)

lower

pore

(

mineral

−

pore

)

−

2(1

−

)(

pore

)

pore

(

pore

)

(2)

where

is the porosity,

pore

and

pore

are the bulk and shear moduli of the pore phase, and

mineral

and

mineral

are the

bulk and shear moduli of the mineral phase.

Figure 3: Schematic of the 3D CNN model. The input to the model is a segmented (pore and mineral) 3D image of size 75

75). The model

is composed mainly of 3D convolution and fully connected linear layers with additional batchnorm, average pooling, non-linear activation layers.

The output of the model is an array of two numbers

and

which lie between 0 and 1. The details of the CNN and fully connected layers are

presented in Table 1

2.3. Details of the training of the surrogate CNN model

The surrogate CNN model is trained on a mixed dataset comprising equal proportion of five sandstone rock images

of size 75

75 voxels. The five sandstone used in this work are the same as those used in Saxena et al. [6]

and Ahmad et al. [8]: two (B1 and B2) are from the Berea Formation, subangular to subrounded Mississippian-age

sandstone; two (FB1 and FB2) from the Fontainebleau Formation, subrounded to rounded Oligocene age sandstone;

and one (CG) from the Castlegate Formation, subangular to subrounded Mesozoic sandstone. All the rocks were

imaged with a micro-CT scanner at the image resolution of approximately 2

m (the physical size of each voxel).

The X-ray di

raction (XRD) analysis demonstrates that all samples mainly consist of quartz mineral along with trace

amounts of feldspar, calcite and clay. For simplicity, we assume single mineralogy where all minerals were treated as

quartz. The rock images are segmented into two phases (pore and mineral) using the Otsu thresholding method [34].

Quartz phase is assigned bulk modulus of 36 GPa and shear modulus of 45 GPa; and pore phase’s bulk and shear

moduli are assumed to be 0 GPa [35] (for numerical stability, pore’s bulk and shear moduli are taken to be 10

−

GPa).

Figure 4 shows the segmented image of rock sample B1 and the two-point correlation function along the

−

direction.

The porosity and correlation function of all five sandstone rocks are presented in Table 2. The rock porosity ranges

from 3.45% to 22.20% and all five rocks have similar correlation length in all three direction indicating their isotropic

structures.

The training data for the surrogate CNN model is generated by homogenizing the subimages by directly solving

the elasticity PDEs using an FFT-based elasticity solver as implemented in GeoDict [36, 37, 38, 39]. The size of the

Table 1: The architecture of the surrogate CNN model. The model comprises two parts: the convolution part and the fully connected part as shown

in Figure 3

The convolution part

Layer

Input size

Kernel size

Number of filters

Output size

Image

Batchnorm3D

Conv3D

PRelu

Dropout3D(prob

0.5)

Batchnorm3D

Conv3D

PRelu

AvgPool3D

Batchnorm3D

Conv3D

PRelu

AvgPool3D

Batchnorm3D

Conv3D

PRelu

Dropout3D(prob

0.5)

Batchnorm3D

Conv3D

PRelu

AvgPool3D

Batchnorm3D

Conv3D

PRelu

Dropout3D(prob

0.5)

The fully connected part

Layer

Input size

Output size

Batchnorm1D

512

Fully connected

PRelu

512

Dropout1D(prob

0.5)

Batchnorm1D

Fully connected

PRelu

Batchnorm1D

Fully connected

0.5(Tanh

Table 2: Statistical properties of the rock samples, including porosity and correlation lengths

, and

directions.

Rock sample

porosity

16.51%

13.1

12.8

13.4

19.63%

14.2

15.3

14.5

22.20%

10.7

8.5

11.8

FB1

9.25%

14.4

16.1

15.6

FB2

3.45%

15.9

15.1

training data is 2800 which contains all five sandstone samples in equal proportion. Moreover, we use 1000 random

rock images from mixed dataset for the validation of the surrogate model. The surrogate CNN model is trained in

supervise manner using the following mean squared loss function

batch

∑

(

∣

DNS

, i

−

pred

, i

∣

DNS

, i

−

pred

, i

∣

)

(3)

B1, porosity = 16.51%

B2, porosity = 19.63%

CG, porosity = 22.20%

FB1, porosity = 9.25%

FB2, porosity = 3.45%

120

160

200

Voxel separation

(

)

0.0

0.2

0.4

0.6

0.8

1.0

Correltation

−

direction

data

fit:

exp

(

−

)

Figure 4: (a) Segmented image of rock sample B1 of size 900

900

900 voxels. Red color represents quartz mineral and white color denotes pore.

(b) Two point correlation function of the segmented image of rock sample B1 along the

−

direction. Blue dots are computed using the segmented

image and the orange curve is the best fit to exponential function exp

(

−

/ξ

)

, with

being the correlation length.

is expressed in the units of

numbers of voxel which is approximately 2

m in length.

where

DNS

and

DNS

are the e

ective bulk and shear moduli obtained from DNS,

pred

and

pred

are the e

ective

bulk and shear moduli predicted by the surrogate CNN model, and

batch

is the batch size. The surrogate CNN model

is trained using the Adam optimizer [40] for 50 epochs where a learning rate of 10

−

is used for the first 25 epochs and

−

for the last 25 epochs. The surrogate CNN is implemented using the PyTorch deep learning framework [41] and

trained on GPU. The variation of the losses during training is shown in Fig. 5(a). The training loss is shown in blue

and the validation loss is shown in orange. Both the training and validation losses decrease during training indicating

that the model is learning to predict e

ective elastic moduli using images. We choose the surrogate CNN model for

which the validation loss s the minimum.

Epoch

0.000

0.025

0.050

0.075

0.100

0.125

0.150

Loss (

)

Training on mixed dataset

Training loss

Validation loss

Epoch

0.0000

0.0025

0.0050

0.0075

0.0100

0.0125

Loss (

)

Transfer learning

FB1

FB2

Figure 5: Variation of the losses (

) during (a) training on mixed dataset and (b) transfer learning. In (a), both training and validation loss functions

are shown. In (b), only the training loss functions are shown for five di

erent sandstones.

After training the surrogate model on the mixed dataset, we make use of

transfer learning

to fine tune the model

for each rock separately. During transfer learning, we fine tune only the fully connected part of the surrogate CNN

model using randomly sampled 1000 images of the rock. The fine-tuning is performed for 50 epochs with a learning

rate of 10

−

for the first 25 epochs and 10

−

for the next 25 epochs. The variation of the training loss functions during

transfer learning are shown in Figure 5(b) for all five rocks. For all five rocks, both the training and validation loss

functions decrease during the transfer learning phase. Thus, after the transfer learning step, we obtain five di

erent

surrogate CNN models, one for each sandstone.

Furthermore, the results obtained from the pretrained CNN model is presented in Appendix A.

3. Results

In this section, we present the results obtained from the application of the CNN-HHM approach(surrogate CNN

model combined with the HHM) to determine e

ective elastic moduli of sandstone images of size 300

300

and 600

600

600 voxels. We also apply the developed framework to limestone with carbonate minerals. All the

results obtained from the CNN-HHM approach are compared with the results obtained from the DNS.