Computing heterogeneous core sample velocity using Digital Rock Physics: A multiscale approach

Journal Pre-proof Computing heterogeneous core sample velocity using Digital Rock Physics: A multiscale approach Sadegh Karimpouli, Asra Faraji, Martin Balcewicz, Erik H. Saenger PII:

S0098-3004(19)30728-9

DOI:

https://doi.org/10.1016/j.cageo.2019.104378

Reference:

CAGEO 104378

To appear in:

Computers and Geosciences

Received Date: 1 August 2019 Revised Date:

15 November 2019

Accepted Date: 18 November 2019

Please cite this article as: Karimpouli, S., Faraji, A., Balcewicz, M., Saenger, E.H., Computing heterogeneous core sample velocity using Digital Rock Physics: A multiscale approach, Computers and Geosciences (2020), doi: https://doi.org/10.1016/j.cageo.2019.104378. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.

1

Computing Heterogeneous Core Sample Velocity Using Using Digital Rock

2

Physics: A Multiscale Approach†

3

Sadegh Karimpoulia*, Asra Farajia, Martin Balcewicz b,c and Erik H. Saenger b,c

4

a. Mining Engineering Group, Faculty of Engineering, University of Zanjan, Zanjan, Iran.

5

* Corresponding author’s Email: [email protected]

6

b. International Geothermal Centre, Bochum University of Applied Sciences, Germany. Email:

7

[email protected]

8

c. Ruhr-University Bochum, Germany. E-mail: [email protected]

9

Abstract

10

Digital Rock Physics (DRP) is an effective approach to compute physical properties of rock using

11

high-resolution 3D images. Although micro-scale structures are well studied in DRP, micro

12

computed tomography (µCT) scanners are still not widely available and imaging could be both

13

expensive and time consuming. Moreover, there is always a trade-off between image resolution

14

and sample size. The later issue is crucial especially in heterogeneous samples such as

15

carbonates. In this study, we propose a multiscale procedure and computed wave velocity of

16

three heterogeneous travertine (calcite) samples with the exact core size of 55 mm diameter.

17

This procedure is conducted using relatively cheap and widely available imaging tools, namely

18

medical CT scanner and conventional microscope. In a low-resolution step, 3D CT-images with

19

200 µm resolution are obtained and used to compute porosity cubes by a three-phase †

Sadegh Karimpouli conceived the problem, developed the idea, performed the coding and contributed to analyzing the results and writing the paper. Asra Faraji prepared core samples, conducted laboratory analysis in Iran and performed image analysis and computation. Martin Balcewicz conducted laboratory analysis in Germany and contributed to writing the paper. Erik Saenger performed computation and contributed to writing the paper.

1

20

segmentation (mineral, pore and sub-resolution). In a high-resolution step, elastic moduli are

21

computed using 2D microscopic images with 1 µm resolution and converted into 3D properties

22

using a 2D-to-3D DRP method. The obtained micro-trends of effective elastic moduli are

23

analytically approximated using the modified Differential Effective Medium (DEM) model. Low-

24

resolution porosity cubes and high-resolution micro-trends are combined to obtain 3D cubes of

25

elastic moduli and density. Finally, a dynamic Finite Difference Method (FDM) is used to

26

propagate P- and S-waves through these samples for velocity computations. Comparison of lab

27

and computed results showed that the errors of P- and S-wave velocities vary from 3.4 % to 7 %

28

and from 6.7 % to 11.1 %, respectively.

29

Keyword: Keyword: Digital Rock Physics (DRP); Heterogeneous carbonate; Multiscale; Core sample velocity

30 31

1. Introduction

32

Knowledge about the wave velocity through the rock is critical tool supporting reservoir

33

evaluation. Velocity and, generally, elastic properties of rock are the main parameters relating

34

seismic data with rock parameters measured in borehole and laboratory. Rock properties are

35

quantified in pore scale using, for example, Digital Rock Physics (DRP) method, which uses high-

36

resolution images of rock and simulates physical processes such as wave propagation to

37

compute velocity. Upscaling of such quantities into well, laboratory (i.e., core) and/or seismic

38

scale is still controversial and, therefore, a direct evaluation of such quantities in a larger scale

39

than pore scale is highly desired. DRP been developed in the last two decades mostly due to

40

emerging of micro X-Ray computed tomography (µCT) scanners and focused ion beam scanning

41

electron microscopy (FIB-SEM) (Andrä et al., 2013a, 2013b). High-resolution representation of

2

42

rock microstructures enabled the DRP to describe the role of complex pore geometry on

43

effective parameters of rock. The standard workflow of DRP is based on three main steps (Andrä

44

et al., 2013a, 2013b; Blunt et al., 2013): (1) imaging, (2) segmentation and (3) numerical

45

computation. The rock sample is imaged using, for example, a µCT scanner which provides a 3D

46

gray-scale image. The reconstructed 3D image is based on a large number of 2D projections in

47

different angles (Mathews et al., 2017). The brightness of image depends on the attenuation

48

coefficient of the material and is strongly linked to density. This means, bright and dark voxels

49

could be representative of mineral and pore phases, respectively. The process of labeling each

50

voxel into different phases is called image segmentation. It is a very critical step in DRP due to its

51

direct effects on pore space structure and, thus, the numerical results. Therefore, several

52

segmentation methods have been introduced including thresholding (Otsu, 1975; Sezgin, 2004),

53

region growing (Beucher and Meyer, 1992; Kass et al., 1988), hybrid (Andrä et al., 2013a;

54

Ramandi et al., 2016) and Deep Learning (DL)-based methods (Karimpouli and Tahmasebi, 2019).

55

Segmentation is even the most critical step when sample fractures and discontinuities are

56

supposed to be delineated (Deb et al., 2008; Karimpouli et al., 2019; Kemeny and Post, 2003;

57

Muller et al., 1992). After image segmentation, proper physical properties could be assigned to

58

each phase and physical processes are numerically computed to simulate elastic wave

59

propagation, electrical current and fluid flow for evaluation of effective elasticity, wave velocity,

60

resistivity and permeability of the sample.

61

Although DRP is an effective method for capturing micro-scale characteristics of rock, there are

62

some limitations. Main limitations are (1) the high-technology µCT scanners are not widely

63

available, (2) imaging is expensive and could be time consuming if higher resolutions are desired,

3

64

(3) sample size is small (mostly in mm order), and (4) a high capacity in terms of computer

65

memory and CPU usage is needed for numerical simulations. The first two issues opened a new

66

avenue through methods, which use 2D this section images to compute 3D properties of rock

67

(Berrezueta et al., 2019; Karimpouli et al., 2018b; Saxena et al., 2017; Saxena and Mavko, 2016).

68

The two last issues are more critical especially in heterogeneous samples such as carbonates,

69

where the sample contains different pores from nano to macroscale. This means, micro porosity

70

could not be resolved by a low-resolution image, but high-resolution image may not be a

71

representative element volume (REV) due to the small field-of-view of the sample (Hebert et al.,

72

2015; Knackstedt et al., 2009). In fact, there is a trade-off between sample size and image

73

resolution. This is one reason why digital models appear to perform stiffer with higher velocities

74

than the measured values of similar samples in laboratory (Arns et al., 2005; Saenger et al.,

75

2011). This difference could also be due to grain-to-grain modulus, crystal deflects and in grain

76

micro-fractures. To overcome this issue and characterize the heterogeneity in different scales,

77

multi-scale imaging is a possible approach (Devarapalli et al., 2017). Although some researchers

78

found similar trends between rock physical parameters in different scales of carbonates (Dvorkin

79

et al., 2011; Saenger et al., 2016), a multi-scale study could be a more reliable approach since it

80

integrates all information from different scales (Hebert et al., 2015). Knackstedt et al. (2009)

81

used thin sections and detected that unresolved pixels in their CT images are due to cemented

82

micrite crystals. They used Self-Consistent (SC) theory to compute elastic moduli of microporous

83

regions and computed more reasonable velocities comparing to experiments. Jouini et al. (2015)

84

proposed a multi-scale approach to extrapolate results from high-resolution (0.3 µm) images

85

into a core scale image with 38 mm diameter. Computed effective moduli by Jouini et al. (2015)

4

86

are overestimated regarding to experimental measurements. They assume nanoscale pores may

87

cause this overestimation, which could be further resolved using nano-tomography. Faisal et al.

88

(2017) tried to directly simulate effective elastic properties of Silurian Dolomite on the core scale

89

(with 38 mm diameter). Their computations in different resolutions showed that effective elastic

90

properties depend not only on the porosity value, but also on microstructure details (we simply

91

call it pore shape). They extracted small non-overlapping cubes from core scale image with

92

different sizes and computed effective elastic moduli of the original sample using elastic moduli

93

of all small cubes. Faisal et al. (2017) showed, that estimations of computed elastic moduli

94

strongly depend on the REVs size, if a constant resolution is set. For instance, a sample with an

95

appropriate REV size results in better estimations than several samples with sizes smaller than

96

the determined REV size. They recently extended their study on similar samples using a multi-

97

scale approach (Farhana Faisal et al., 2019), where images with three resolutions were used: 13

98

µm (low-resolution), 1 µm (medium-resolution), and 0.12 µm (high-resolution). To resolve the

99

sub-resolution parts of the low-resolution image, they prepared limited numbers of images with

100

medium-resolution. They also used some high-resolution FIB-SEM images to resolve nanoscale

101

pores of the medium images. Effective moduli of small images were computed, and their

102

average values were used as the moduli of the sub-resolution parts of the medium-resolution

103

image. Similarly, the average properties of medium-resolution images were used for sub-

104

resolution parts of the low-resolution image. Results showed reasonable improvements, where

105

overestimated errors were reduced from 8.9 to 4.5 % for P-wave velocity and from 11.9 to 7.8 %

106

for S-wave velocity. Although, these works represented valuable results, they neglected the

107

effect of pore shapes. In fact, in two latter researches, Faisal et al. (2017, 2019) determined that

5

108

microstructural details like pore shape and structure play an important role, but they suppressed

109

this effect by choosing an average value. We suggest instead of using the average properties for

110

the unresolved part, microscale trends of elastic moduli should be studied and used in

111

macroscale images. As a similar work, Karimpouli et al. (2018a) used such an idea to resolve the

112

sub-resolution carbonate cement in a cemented sandstone. Due to limitation of the available

113

data, they used a random procedure for porosity value of sub-resolution pixels. Subsequently,

114

they found two negative effects: (1) the spatial distribution of sub-resolution porosity is

115

neglected and (2) the overall properties of sub-resolution pixels approaches to average

116

properties of the medium (similar to Faisal et al. (2017, 2019)). This is why elastic moduli were

117

also overestimated in their study. Unlike to that work, we here demonstrate, this idea could lead

118

to better results if a reliable procedure is used for both porosity evaluation of the sub-resolution

119

pixels and assigning proper elastic properties based on elastic property-porosity trends in

120

microscale.

121

High resolution imaging requires small subsamples, if µCT or FIB-SEM methods are applied.

122

Preparing subsamples results into numerous recordings which are expensive, time consuming,

123

and not widely available. Therefore, we aim to directly compute velocity of a core sample (e.g.,

124

with about 5 cm diameter and 10 cm length) and tackle the multiscale problem using

125

inexpensive and widely spread imaging tools. Our methodology is subdivided into three steps.

126

First, we use a medical CT scanner, which is both cheap and accessible for the macro-imaging.

127

The limitation for this technology is pore size. The size of the largest pores in the sample must be

128

at least two-times more than the CT resolution. In a second step, we use 2D microscopic image

129

and try to convert 2D results into 3D characteristics using alternative DRP methods (Saxena and

6

130

Mavko, 2016). Finally, we merge high resolution 3D characteristics with lower resolution medical

131

CT images to calculate effective properties of the rock.

132

The rest of the paper is organized as follows. In section 2, core samples and their lab

133

measurements are introduced. Then, the multiscale procedure is explained in section 3. Section

134

4 presents computations and results and, finally, section 5 summarizes the limitations and

135

advantages of the methodology.

136 137

2. Core samples

138

In this study, we use a travertine (limestone) sample from west of Zanjan, Iran. Travertine layers

139

were formed by calcareous hot springs originated from tertiary volcanic intrusions. It is a highly

140

heterogeneous sample with pores sizes varying from cm to even nm. This lithology has been

141

chosen due to macro- and microscopic heterogeneity, and to demonstrate our multiscale

142

procedure using medical CT scanners. Three core samples, namely S1, S2, and S3, were drilled

143

with 5.5 cm diameter and 10 cm length, cut plane-parallel and polished (Fig. 1). After

144

preparation all samples were oven-dried at 60 °C for 48 hours. As it is observed in this figure,

145

each sample contains large (about 1-2 cm) and micro pores (less than 1 µm). Spatial distribution

146

of pore spaces is also variable from well distributed (S2) to poorly distributed (S3) sample.

147

Chemical analysis (Table 1) showed that more than 99 % composition of travertine sample is

148

calcite. Physical characteristics of these samples were measured in the laboratory and

149

summarized in Table 2. Porosimetry was conducted using helium injection method and wave

150

velocity was measured by ultrasonic pulse velocity test at ambient pressure-temperature

7

151

conditions in the laboratory of Research Institute of Petroleum of Industry (RIPI), Iran. Wave

152

velocity was also measured in laboratory of International Geothermal Center (GZB), Germany.

153

Ultrasonic measurements were conducted by similar set-ups at RIPI and GZB. For example, in

154

Germany, a table-top solution was used. The GZB set-up consists mainly of three parts. A

155

waveform generator, two identical ultrasonic sensors (transmitter/receiver) and an oscilloscope.

156

The external waveform generator (Panametrix EPOCH 650) creates a rectangular signal (400 V)

157

with a controlled central frequency (1 MHz), which causes the piezoelectric transducer to emit

158

mechanical pulses into the sample. At the other sample's end, the receiver converts the arriving

159

mechanical pulses (waveforms) back into electric signals. The digital oscilloscope (Picoscope

160

5444B) records the electric waveforms. P- and S-wave arrival signals were determined by manual

161

picking. Effective travel times for each sample was ensured by correcting travel times in the set-

162

up. P- and S-wave velocities were detected for dry state.

163 164

3. Multiscale procedure

165

As it is illustrated in Fig. 2, the heterogeneous travertine sample contains pore spaces with a size

166

from cm to µm or even nm. This means, in a low-resolution image, macro-pores are detected,

167

and sub-resolution micro-pores are missed. In a similar manner, a high-resolution image resolves

168

micro-pores, but macro-pores are missed due to small field-of-view.

169

Our suggested multiscale procedure is visually explained in Fig. 3. In this procedure, we use two

170

high- and low-resolution images, but one can extend this procedure with more images in

171

different resolutions. Large size pores of travertine made it reasonable for us to use medical CT

172

scanners for capturing the low-resolution images with all relevant features. The largest cube

8

173

from the center of the core is extracted as the digital image of each sample. The size of this cube

174

is relatively small (e.g. less than 5003 voxels) even for a sample in cm scale. This is a significant

175

matter from a computational standpoint, where numerical computations could be also

176

conducted very cheap. In this step, a three-phase segmentation (Jouini et al., 2015) is used to

177

divide the CT-image into a porosity cube. In this cube, porosity of each voxel is either zero, if it is

178

a mineral phase, 100 %, if it is a pore phase, or a value between 0 to 100 % depending on gray-

179

scale level of the sub-resolution voxel. Then, material-specific properties (i.e. elastic moduli and

180

density) are assigned to each phase. For the mineral and pore phases with 0 and 100 %,

181

corresponding mineral and fluid properties (e.g. air in dry condition) were applied. However, the

182

question is: “Which values should be assigned to sub-resolution voxels?” A widespread method is

183

to use average values computed from high-resolution images with micro-scale pores (Faisal et

184

al., 2017; Farhana Faisal et al., 2019; Jouini et al., 2015), which suppresses the effects of micro-

185

scale pore shape and value. Previous studies in carbonates showed, that the effects of pore

186

shape could be as important as pore values for evaluation of velocity (Anselmetti and Eberli,

187

1993; Eberli et al., 2003; Karimpouli et al., 2013). This means, effective properties of sub-

188

resolution voxels should be obtained using real trends, which are behavior of one physical

189

property versus another property of rock (e.g. bulk modulus versus porosity). Since we obtain

190

these trends in micro-scale, we call them “micro-trends”, which are affected by exact micro-pore

191

values and structures.

192

To find micro-trends, high-resolution images are needed. These images could be captured either

193

by µCT scanners or FIB-SEM, but for the sake of a simple and cost-effective procedure, we use

194

2D microscopic images. High-resolution images are captured from polished thin sections (Fig. 3).

9

195

The main point is, that the size of high-resolution images should be equal to one low-resolution

196

voxel length. Generally, it is not necessary that the size of high-resolution image is exactly equal

197

to the voxel-size of the low resolution image, but the best state is theoretically when that are

198

equal.

199

In this case, these images are segmented into two phases, namely mineral and pore. After

200

segmenting, 2D effective elastic properties could be computed by either finite element (FEM) or

201

finite difference methods (FDM) codes. The solver related uncertainty is in order of 1/10th of

202

input mineral moduli (Saxena et al., 2019) and could be ignored compared with other major

203

controls such as image resolution. To convert 2D into 3D effective properties, several different

204

DRP methods could be used (Karimpouli et al. 2018b; Karimpouli and Tahmasebi 2016; Saxena

205

and Mavko 2016). Analytical effective models such as the Weighted Average Boundary (WAB)

206

model or the Differential Effective Medium (DEM) model describe micro-trends. Since these

207

models are fitted onto effective moduli computed by exact micro-scale structures of the sample,

208

they enact the effects of pore shape and value simultaneously.

209

Micro-trends can be used to assign to sub-resolution voxels of the low-resolution porosity cube

210

appropriate effective elastic properties (Fig. 3). This leads to cubes of elastic moduli and density

211

from which the effective velocities are computed. In practice, FEM code (Garboczi and Day,

212

1995) can be applied if there are less than 100 phases in the image (Garboczi, 1998), however, in

213

this procedure each voxel is considered as one phase due to its individual properties. This means,

214

for a 3D cube with voxel numbers in order of millions, FEM code is impractical. The used dynamic

215

FDM code (Saenger et al., 2000; Saenger and Shapiro, 2002) does not struggle with allocation

216

limitations.

10

217

3.1. Effective elastic models

218

Estimation of effective elastic moduli of a mixture of grains and pores deals with the volume

219

fraction, elastic moduli, and geometrical details of each phase. By specifying volume fraction and

220

elastic moduli, we can predict upper (stiffer) and lower (softer) material bounds (Mavko et al.,

221

2009). Hashin and Shtrikman (1963) represented the best bounds for an isotropic linear elastic

222

composite known as Hashin-Shtrikman bounds. Dvorkin and Nur (1996) introduced stiff and soft

223

sand models for cemented and uncemented sands, but were also applied for general trends of

224

carbonates (Andrä et al., 2013b). Bounding models reveal the effective elastic properties of two

225

endmembers, but for a sample located between two bounds, Weighted Average Bounding

226

(WAB) model could be used. To consider the third parameter, geometrical details, the methods

227

based on Self-consistent (SC) or Differential Effective Medium (DEM) theory could be used.

228

These models are based on pore shape and generate reliable effective elastic moduli according

229

to their fields of applications (Wang et al., 2015; Xu and Payne, 2009).

230

3.1.1. Differential Effective Medium (DEM)

231

In DEM theory, two-phase composites are modeled by considering one host (or mineral) phase

232

and incrementally adding guest (or pore) phase. When porosity is zero, mineral is assumed as

233

the matrix or host phase. Porosity with known geometrical details or pore shape is added

234

incrementally until a target porosity value is obtained. In each step, effective elastic moduli are

235

computed and used as the matrix elastic moduli for the next step (Mavko et al., 2009). The DEM

236

equations are:
1 −   ∗
 =
− ∗ 
∗ 
,  

(1)

11


1 −  237 238 239





 ∗
 =
 −  ∗ 
∗ 
,

(2)

where ∗ and  ∗ are effective elastic moduli with initial conditions ∗ =  and  ∗ =  ,

where  and  are the bulk and shear moduli of the host material (mineral), and  are the

bulk and shear moduli of the incrementally added inclusions,  is the concentration of guest

240

phase (porosity) and 
∗  and 
∗  are geometrical factors of guest phase (for more details see

241

Mavko et al. (2009)). Geometrical factors determine the effects of guest inclusions with different

242

materials and shapes. The computational details are explained in Appendix A. For a two-phase

243

medium, as in our case, geometrical factors are obtained according to only the shape of

244

inclusions. Each inclusion is assumed as an ellipsoid and the corresponding ratio of small to large

245

radii, known as aspect ratio, is used as the shape factor to compute geometrical factors. Effective

246

elastic moduli are highly affected by the aspect ratio of pore spaces (Berryman, 1980). In a

247

randomly oriented model, it is stiff model if the aspect ratio is one (spherical pores), and it

248

becomes softer if the aspect ratio is reduced (ellipsoidal pores). Where the aspect ratio is very

249

small, pore space is like a fracture and a weak model is obtained. In DEM model, porosity can

250

increase up to 100 % if solid phase is connected. However, in rock samples, there is a critical

251

porosity ( ), where solid matrix falls apart. Mukerji et al. (1995) proposed modified DEM and

252

assumed percolation behavior of the rock at critical porosity as the guest phase. Therefore, in

253

modified DEM,  denotes the concentration of the critical phase in the matrix.

254 255

4. Results

256

4.1. Low-resolution step

12

257

All core samples are imaged by a medical CT scanner, which are shown in Fig. 4. The largest cube

258

is extracted from each image called “CT-image” (Fig. 4). For each sample, the size of the CT

259

image is 200×200×475 voxels, where edge length of each voxel is 200 µm. According to image

260

resolution all pores with a size of less than 200 µm are not resolved and the containing voxels

261

represent a gray-scale value depending on the amount of occupied grains and pores.

262

To have an insight about porosity and pore structure of the samples, large-scale pores were

263

divided using a two-phase segmentation by selecting a rough threshold value. Fig. 4 illustrates

264

the pore structure and distribution through 3D space of the samples. First hand computations

265

(i.e., two phase segmentation) showed that porosity of samples is roughly about 15 % for

266

samples S1 and S2, and 10 % for S3. These values are not consistent with laboratory

267

measurements (Table 2), referring 9.8 % for S1 and S2, and 7.4 % for S3, respectively. Fig. 4

268

reveals that there are many large and isolated pores with no connection with connected pore

269

system, which are not measured during helium-injection porosimetry. According to the nature of

270

travertines, when hot spring waters rich in calcium carbonate and dissolved carbon dioxide are

271

exposed to Earth’s surface, carbon dioxide is outgassed and calcite (and occasionally aragonite)

272

may precipitate from supersaturated solution of carbonate calcium (Drysdale, 1999). Gasification

273

could lead to isolated pores with any size from cm to nm. Distribution of such pores could also

274

be very heterogeneous due to variable precipitation process of travertines.

275

To explore if the CT-image is larger than the Representative Volume Element (RVE), in each

276

sample, several cubes with different sizes were extracted from the center of CT-image and

277

corresponding porosity were computed. Porosity variation with image size is plotted in Fig. 5.

13

278

According to this plot, almost in all sample, a cube with a size of 1503 voxels could be considered

279

as RVE.

280

Since the resolution of CT-images is low, two-phase segmentation could not produce reliable

281

results. For tacking sub-resolution porosity into account, three-phase segmentation is used (Fig.

282

6) (Jouini et al., 2015). This method uses two threshold values ( and  in Fig. 6) and

283

divides all voxels into three phases: mineral, pore and sub-resolution phases. The corresponding

284

porosity of mineral and pore phases are 0 and 100 %, but the porosity of sub-resolution phase is

285

computed by a simple linear interpolation between two threshold values (Fig. 6). However, a

286

more reliable way for computing porosity of sub-resolution phase is to use a high-resolution

287

image. Supposing, if all micro-structures are resolved in such an image, an empirical relation is

288

obtained by correlating micro-porosity and gray-scale values (Ramandi et al., 2016). In this study,

289

we used the linear interpolation. One can use K-means clustering (MacQueen, 1967) or bi-level

290

segmentation (Ridler and Calvard, 1978) methods for finding threshold values, but Jouini et al.

291

(2015) set them manually for more accurate results. We also used a porosity inverse method

292

(Yao et al., 2009) to set the threshold values. To this end, porosity of each sample is obtained

293

using bulk and grain density measurements, and threshold values are set accordingly. Porosity

294

values were computed as 15.04 %, 14.29 % and 10.95 % for S1, S2, and S3, respectively. The

295

porosity cubes could be used to assign corresponding values of elastic moduli and density to

296

each phase. We used properties of calcite and air for mineral and pore phases (Table 3), but for

297

the sub-resolution phase micro-trends are needed, which are obtained using high-resolution

298

scale.

299

4.2. High-resolution step

14

300

In this step, four polished thin sections were prepared and imaged using conventional optical

301

microscopes. These sections were not select from the cores, which means our high-resolution

302

images are not registered to low-resolution ones. The size of high-resolution images is set to 200

303

µm with a resolution of 1 µm/pixel, which is equal to low-resolution voxel length (200 µm/voxel).

304

Since 2D imaging and corresponding numerical processing are fast, 200 individual 2D high-

305

resolution microscopic images were captured (Fig. 7). The microscope settings and

306

corresponding imaging software were arranged to show high contrast. Reflective light was

307

applied to indicate pores as black spots and minerals as grey areas. This makes it easy to

308

segment them using a simple thresholding. We used Otsu thresholding (Otsu, 1975) to segment

309

the grain and pore phases. Results of segmentation are shown in Fig. 7.

310

To compute effective elastic properties from 2D images of rock samples, the Finite Element code

311

by Garboczi and Day (1995) was used, which is based on solving the Hooke’s law equations of

312

linear elasticity. Results are shown in Fig. 8. These results are in 2D and, before any further

313

process, they should be converted into 3D properties. By using 2D information to predict 3D rock

314

properties we assume a similar heterogeneity level in all three spatial directions. Among all

315

alternative DRP methods, we used the method by Saxena and Mavko (2016), which converts

316

elastic properties from 2D into 3D using the following equations: 

=

"

=



"

317 318

 ,

(3)

 ,

(4)

 

" "

!

#

where and  are bulk and shear moduli of three- and two-dimensional models (3% and 2%) and mineral ('). ( and (" are constant powers, which should be adjusted in any case. Saxena

15

319 320

and Mavko (2016) proposed ( = (" = 0.6 for carbonate samples. Also, we found ( =

0.6, (" = 0.5 for real carbonate reservoir rocks (Karimpouli et al., 2018a). In this case, we used

321

an average value of core measured velocity and determined ( = (" = 0.7 are reasonable

322

values. Results of 3D effective elastic moduli of high-resolution images are shown in Fig. 7. These

323

results reveal the microscopic properties of the sample, which are affected by real pore shapes

324

and structures.

325

To find micro-trends, effective elastic model based on DEM equations (Eq. 1, 2) were used. In the

326

DEM model, bulk and shear aspect ratios (needed for 
∗  and 
∗  in Eq. 1, 2) should be

327

optimized for obtaining the best fit. The criterion for best fitting is Mean Square Error (/0): /0 =

328



∑= 4>?234567 839:7; < 

(5)

,

where / @ and /A@B are effective elastic moduli computed from images (we mean the 3D

329

values) and calculated by models and C is the number of images. In modified DEM, two aspect

330

ratios are obtained for bulk and elastic moduli. The optimum aspect ratios were obtained as 0.74

331

and 0.21 with MSE values of 1072 and 1982 (GPa)2 for modified DEM models of bulk and shear

332

moduli. These micro-trends are shown in Fig. 7. In this figure, porosity axes are plotted up to 60

333

% porosity, however, the micro-trends extend up to 100 % porosity where they are 0.

334

4.3. Combined high- and low-resolution step

335

The micro-trends (DEM models) obtained in the last step could be used to assign proper

336

properties to each voxel of low-resolution image based on the porosity cube. For each voxel, if

337

porosity is 0 or 100 %, properties of calcite or air (Table 3) are assigned, respectively. For other

338

porosity values, properties are obtained from the modified DEM trends and assigned into the

16

339

corresponding voxels, which leads to low-resolution cubes of elastic moduli. Since porosity of

340

each voxel is known, the following relation is used to compute cube of density: DEFBG = DH +
1 − D ,

(6)

341

where DH and D are density of fluid (here air) and mineral.

342

Finally, the dynamic pulse propagation method (e.g. Saenger and Shapiro (2002)), which is based

343

on a staggered-grid FDM, is used to compute P- and S-wave velocities using cubes of elastic and

344

density properties. The heterogeneous sample constructed by these cubes is embedded in a

345

homogeneous elastic region with elastic properties of calcite. A plane-source wave is generated

346

and propagated numerically through the sample using a periodic boundary condition. The time-

347

delay of the peak amplitude of the mean plane wave could be used to estimate effective P- and

348

S-wave velocities. Results for velocity computations are plotted in Fig. 9 and summarized in Table

349

4.

350 351

5. Discussions

352

According to our results (Table 4), the mean percentage errors (Δ %) of P-wave velocity, defined

353

as the difference between lab and computed velocity, are 3.4 % and 7 % for RIPI and GZB

354

laboratories. The errors of S-wave velocity are 11.1 % and 6.7 % for RIPI and GZB laboratories.

355

The smallest errors of the alternative multiscale studies (Faisal et al., 2017; Farhana Faisal et al.,

356

2019; Jouini et al., 2015) are 4.5 and 7.8 % for P- and S-wave velocity. Comparing to these

357

studies, which are expensive in both imaging tools and computation points of view, our results

358

are reasonable and solid.

17

359

There are several parameters, which impact the results. For example, medical CT is a proper tool

360

for porous carbonates such as travertine, but it could not be used for other reservoir rocks with

361

pore sizes smaller than about 200 µm, for instance, most of sandstones and compacted

362

carbonates. However, the proposed multiscale procedure (Fig. 3) could be used by other

363

techniques on different scales. It means, instead of using images in mm and µm scales such as

364

what we had in our case, one can use some images in µm and nm scales with similar multiscale

365

procedure. This procedure could also be extended in more than two steps. The main points for

366

such a procedure are: (1) the image size of each step should be equal to the voxel size of the last

367

step, and (2) segmentation of all steps are three-phase except for the highest-resolution step,

368

which is two-phase segmentation.

369

We limited ourselves to use 2D microscopic images in the high-resolution step, because we

370

wanted to find the most economical but still accurate way. Using 3D images in this step could

371

lead to more reliable results depending on rock type and heterogeneity. In our study, we use a

372

huge number of images (i.e., 200) for the case of carbonates. The resulting rock physical trend

373

should be representative even for this relative heterogeneous rock type. In addition, these

374

images were not registered. The authors suppose that being registered does not have a drastic

375

effect on final results, however, we suggest high-resolution images to be registered if smaller

376

errors are desired.

377

Finally, the key point of this multiscale study is using micro-trends instead of an average value for

378

sub-resolution pixels. The effective properties model, used to fit the micro-trends, could be

379

variant in different samples. Since they are computed from the real micro-scale world of the

380

sample, they are affected by the of pore shapes of the sample. It means, the effects of pore

18

381

shape, which is considered as one of the most important rock physics parameters in carbonates,

382

are implied in numerical computation, which leads to more realistic results.

383 384

6. Conclusion

385

In this study, three core samples with 5.5 cm diameter and 10 cm length from a highly

386

heterogeneous travertine were selected for velocity computation using DRP. We introduced a

387

multiscale procedure using relatively cheap and widely available imaging tools. In the low-

388

resolution step, we used a medical CT tool for 3D imaging, which generates images with the size

389

reasonable for a fast-numerical computation. Since micro-pores are not resolved in this scale,

390

three-phase segmentation was used to divide CT-images into mineral, pore and sub-resolution

391

voxels. To assign physical properties to each voxel, elastic properties of calcite and air properties

392

were used for voxels with mineral and pore phases. However, for sub-resolution voxels, we used

393

two threshold values and predicted porosity using a linear interpolation. Corresponding

394

properties to each porosity were obtained from micro-trends of the sample. In the high-

395

resolution step, conventional optical microscopes were applied to capture 2D images in order to

396

show the ability of alternative DRP methods. Micro-trends were modeled by WAB and modified

397

DEM models, but DEM was better fitted with smaller MSE for both bulk and shear moduli. Since

398

these trends capture micro-scale pore shape and structure, they highly improve the estimation

399

of effective elastic moduli. Using these micro-trends, cubes of effective elastic moduli and

400

density were computed based on low-resolution CT-image and porosity cube. Finally, a FDM

401

code was applied to compute effective wave velocity. Our results showed that the errors of P-

19

402

and S-wave velocity computations are 3.4 % and 11.1 % for RIPI measurements and they are 7 %

403

and 6.7 % for GZB measurements.

404

The proposed procedure was used for a sample with large pores. Though, it could be

405

implemented for samples with pore sizes smaller than the resolution of medical CT scanners. In

406

this case, images with smaller scales are needed. This procedure could also be extended to more

407

than two scale levels. For a routine use of such a multi-scale approach some basic calibration

408

steps are necessary as we discussed in this paper.

409 410

References

411

Andrä, H., Combaret, N., Dvorkin, J., Glatt, E., Han, J., Kabel, M., Keehm, Y., Krzikalla, F., Lee, M.,

412

Madonna, C., Marsh, M., Mukerji, T., Saenger, E.H., Sain, R., Saxena, N., Ricker, S.,

413

Wiegmann, A., Zhan, X., 2013a. Digital rock physics benchmarks-Part I: Imaging and

414

segmentation. Comput. Geosci. 50, 25–32. https://doi.org/10.1016/j.cageo.2012.09.005

415

Andrä, H., Combaret, N., Dvorkin, J., Glatt, E., Han, J., Kabel, M., Keehm, Y., Krzikalla, F., Lee, M.,

416

Madonna, C., Marsh, M., Mukerji, T., Saenger, E.H., Sain, R., Saxena, N., Ricker, S.,

417

Wiegmann, A., Zhan, X., 2013b. Digital rock physics benchmarks-part II: Computing effective

418

properties. Comput. Geosci. 50, 33–43. https://doi.org/10.1016/j.cageo.2012.09.008

419

Anselmetti, F.S., Eberli, G.P., 1993. Controls on sonic velocity in carbonates. Pure Appl. Geophys.

420

PAGEOPH 141, 287–323. https://doi.org/10.1007/BF00998333

421

Arns, C.H., Bauget, F., Limaye, A., Sakellariou, A., Senden, T., Sheppard, A., Sok, R.M., Pinczewski,

422

V., Bakke, S., Berge, L.I., Oren, P., Knackstedt, M., 2005. Pore-Scale Characterization of

423

Carbonates Using X-Ray Microtomography. SPE J. 10, 26–29.

20

424 425

https://doi.org/10.2118/90368-PA Berrezueta, E., Domínguez-Cuesta, M.J., Rodríguez-Rey, Á., 2019. Semi-automated procedure of

426

digitalization and study of rock thin section porosity applying optical image analysis tools.

427

Comput. Geosci. 124, 14–26. https://doi.org/10.1016/j.cageo.2018.12.009

428 429 430 431 432 433 434

Berryman, J.G., 1980. Long-wavelength propagation in composite elastic media II. Ellipsoidal inclusions. J. Acoust. Soc. Am. 68, 1820–1831. Beucher, S., Meyer, F., 1992. The morphological approach to segmentation: the watershed transformation. Opt. Eng. YORK-MARCEL DEKKER Inc. 34, 433–433. Blunt, M.J., Bijeljic, B., Dong, H., Gharbi, O., Iglauer, S., Mostaghimi, P., Paluszny, A., Pentland, C., 2013. Pore-scale imaging and modelling. Adv. Water Resour. 51, 197–216. Deb, D., Hariharan, S., Rao, U.M., Ryu, C.H., 2008. Automatic detection and analysis of

435

discontinuity geometry of rock mass from digital images. Comput. Geosci. 34, 115–126.

436

https://doi.org/10.1016/j.cageo.2007.03.007

437

Devarapalli, R.S., Islam, A., Faisal, T.F., Sassi, M., Jouiad, M., 2017. Micro-CT and FIB–SEM imaging

438

and pore structure characterization of dolomite rock at multiple scales. Arab. J. Geosci. 10,

439

361. https://doi.org/10.1007/s12517-017-3120-z

440

Drysdale, R.N., 1999. The sedimentological significance of hydropsychid caddis-fly larvae (order;

441

Trichoptera) in a travertine-depositing stream; Louie Creek, Northwest Queensland,

442

Australia. J. Sediment. Res. 69, 145–150. https://doi.org/10.2110/jsr.69.145

443 444 445

Dvorkin, J., Derzhi, N., Diaz, E., Fang, Q., 2011. Relevance of computational rock physics. Geophysics 76, E141–E153. https://doi.org/10.1190/geo2010-0352.1 Dvorkin, J., Nur, A., 1996. Elasticity of high-porosity sandstones: Theory for two North Sea data

21

446 447

sets. GEOPHYSICS 61, 1363–1370. https://doi.org/10.1190/1.1444059 Eberli, G.P., Baechle, G.T., Anselmetti, F.S., Incze, M.L., 2003. Factors controlling elastic

448

properties in carbonate sediments and rocks. Lead. Edge 22, 654–660.

449

https://doi.org/10.1190/1.1599691

450

Faisal, T.F., Awedalkarim, A., Chevalier, S., Jouini, M.S., Sassi, M., 2017. Direct scale comparison

451

of numerical linear elastic moduli with acoustic experiments for carbonate rock X-ray CT

452

scanned at multi-resolutions. J. Pet. Sci. Eng. 152, 653–663.

453

https://doi.org/10.1016/j.petrol.2017.01.025

454

Farhana Faisal, T., Islam, A., Jouini, M.S., Devarapalli, R.S., Jouiad, M., Sassi, M., Faisal, T.F., Islam,

455

A., Jouini, M.S., Devarapalli, R.S., Jouiad, M., Sassi, M., 2019. Numerical prediction of

456

carbonate elastic properties based on multi-scale imaging. Geomech. Energy Environ.

457

100125.

458

Garboczi, E.J., 1998. Finite Element and Finite Difference Programs for Computing the Linear

459

Electric and Elastic Properties of Digital Images of Random Materials | NIST. NIST

460

Interagency/Internal Rep. - 6269.

461

Garboczi, E.J., Day, A.R., 1995. An algorithm for computing the effective linear elastic properties

462

of heterogeneous materials: three-dimensional results for composites with equal phase

463

Poisson ratios. J. Mech. Phys. Solids 43, 1349–1362. https://doi.org/10.1016/0022-

464

5096(95)00050-S

465 466 467

Hashin, Z., Shtrikman, S., 1963. A variational approach to the theory of the elastic behaviour of multiphase materials. J. Mech. Phys. Solids 11, 127–140. Hebert, V., Garing, C., Luquot, L., Pezard, P.A., Gouze, P., 2015. Multi-scale X-ray tomography

22

468

analysis of carbonate porosity. Geol. Soc. London, Spec. Publ. 406, 61–79.

469

https://doi.org/10.1144/SP406.12

470

Jouini, M.S., Vega, S., Al-Ratrout, A., Al-Ratrout, A., 2015. Numerical estimation of carbonate rock

471

properties using multiscale images. Geophys. Prospect. 63, 405–421.

472

https://doi.org/10.1111/1365-2478.12156

473

Karimpouli, S., Hassani, H., Nabi-Bidhendi, M., Khoshdel, H., Malehmir, A., 2013. Application of

474

probabilistic facies prediction and estimation of rock physics parameters in a carbonate

475

reservoir from Iran. J. Geophys. Eng. 10, 1–14. https://doi.org/10.1088/1742-

476

2132/10/1/015008

477

Karimpouli, S., Khoshlesan, S., Saenger, E.H., Koochi, H.H., 2018a. Application of alternative

478

digital rock physics methods in a real case study: a challenge between clean and cemented

479

samples. Geophys. Prospect. Accepted, Published Online. https://doi.org/10.1111/1365-

480

2478.12611

481 482 483 484 485

Karimpouli, S., Tahmasebi, P., 2019. Segmentation of Digital Rock Images Using Deep Convolutional Autoencoder Networks. Comput. Geosci. 126, 142–150. Karimpouli, S., Tahmasebi, P., 2016. Conditional reconstruction: An alternative strategy in digital rock physics. GEOPHYSICS 81, D465–D477. https://doi.org/10.1190/geo2015-0260.1 Karimpouli, S., Tahmasebi, P., Saenger, E.H., 2019. Coal Cleat/Fracture Segmentation Using

486

Convolutional Neural Networks. Nat. Resour. Res. https://doi.org/10.1007/s11053-019-

487

09536-y

488 489

Karimpouli, S., Tahmasebi, P., Saenger, E.H., 2018b. Estimating 3D elastic moduli of rock from 2D thin-section images using differential effective medium theory. GEOPHYSICS 83, MR211–

23

490 491 492 493

MR219. https://doi.org/10.1190/geo2017-0504.1 Kass, M., Witkin, A., Terzopoulos, D., 1988. Snakes: Active contour models. Int. J. Comput. Vis. 1, 321–331. https://doi.org/10.1007/BF00133570 Kemeny, J., Post, R., 2003. Estimating three-dimensional rock discontinuity orientation from

494

digital images of fracture traces. Comput. Geosci. 29, 65–77.

495

https://doi.org/10.1016/S0098-3004(02)00106-1

496

Knackstedt, M.A., Latham, S., Madadi, M., Sheppard, A., Varslot, T., Arns, C., 2009. Digital rock

497

physics: 3D imaging of core material and correlations to acoustic and flow properties. Lead.

498

Edge 28, 28–33. https://doi.org/10.1190/1.3064143

499

MacQueen, J.B., 1967. Some Methods for classification and Analysis of Multivariate

500

Observations, in: Proceedings of 5th Berkeley Symposium on Mathematical Statistics and

501

Probability. University of California Press, pp. 281–297.

502

Mathews, J.P., Campbell, Q.P., Xu, H., Halleck, P., 2017. A review of the application of X-ray

503

computed tomography to the study of coal. Fuel 209, 10–24.

504

https://doi.org/10.1016/J.FUEL.2017.07.079

505

Mavko, G., Mukerji, T., Dvorkin, J., 2009. The rock physics handbook: Tools for seismic analysis of

506

porous media. Cambridge university press.

507

https://doi.org/http://dx.doi.org/10.1017/CBO9780511626753

508

Mukerji, T., Berryman, J., Mavko, G., Berge, P., 1995. Differential effective medium modeling of

509

rock elastic moduli with critical porosity constraints. Geophys. Res. Lett. 22, 555–558.

510

https://doi.org/10.1029/95GL00164

511

Muller, J.J., Schrauber, H., Müller, J.J., Schrauber, H., 1992. The inertia-equivalent ellipsoid: a link

24

512

between atomic structure and low-resolution models of small globular proteins determined

513

by small-angle X-ray scattering. J. Appl. Crystallogr. 25, 181–191.

514

https://doi.org/10.1107/S0021889891011421

515

Otsu, N., 1975. A threshold selection method from gray-level histograms. Automatica 11, 23–27.

516

Ramandi, H.L., Armstrong, R.T., Mostaghimi, P., 2016. Micro-CT image calibration to improve

517

fracture aperture measurement. Case Stud. Nondestruct. Test. Eval.

518

Ridler, T.W., Calvard, S., 1978. Picture Thresholding Using an Iterative Selection Method. IEEE

519

Trans. Syst. Man. Cybern. 8, 630–632. https://doi.org/10.1109/TSMC.1978.4310039

520

Saenger, E.H., Enzmann, F., Keehm, Y., Steeb, H., 2011. Digital rock physics: Effect of fluid

521

viscosity on effective elastic properties. J. Appl. Geophys. 74, 236–241.

522

https://doi.org/10.1016/j.jappgeo.2011.06.001

523

Saenger, E.H., Gold, N., Shapiro, S.A., 2000. Modeling the propagation of elastic waves using a

524

modified finite-difference grid. Wave Motion 31, 77–92. https://doi.org/10.1016/S0165-

525

2125(99)00023-2

526

Saenger, E.H., Shapiro, S.A., 2002. Effective velocities in fractured media: a numerical study using

527

the rotated staggered finite-difference grid. Geophys. Prospect. 50, 183–194.

528

https://doi.org/10.1046/j.1365-2478.2002.00309.x

529

Saenger, E.H., Vialle, S., Lebedev, M., Uribe, D., Osorno, M., Duda, M., Steeb, H., 2016. Digital

530

carbonate rock physics. Solid Earth 7, 1185–1197. https://doi.org/10.5194/se-7-1185-2016

531

Saxena, N., Hofmann, R., Hows, A., Saenger, E.H., Duranti, L., Stefani, J., Wiegmann, A., Kerimov,

532

A., Kabel, M., 2019. Rock compressibility from microcomputed tomography images:

533

Controls on digital rock simulations. GEOPHYSICS 84, WA127–WA139.

25

534 535

https://doi.org/10.1190/geo2018-0499.1 Saxena, N., Mavko, G., 2016. Estimating elastic moduli of rocks from thin sections: Digital rock

536

study of 3D properties from 2D images. Comput. Geosci. 88, 9–21.

537

https://doi.org/10.1016/j.cageo.2015.12.008

538

Saxena, N., Mavko, G., Hofmann, R., Srisutthiyakorn, N., 2017. Estimating permeability from thin

539

sections without reconstruction: Digital rock study of 3D properties from 2D images.

540

Comput. Geosci. 102, 79–99. https://doi.org/10.1016/j.cageo.2017.02.014

541 542 543

Sezgin, M., 2004. Survey over image thresholding techniques and quantitative performance evaluation. J. Electron. Imaging 13, 146–168. Wang, Z., Wang, R., Li, T., Qiu, H., Wang, F., 2015. Pore-scale modeling of pore structure effects

544

on P-wave scattering attenuation in dry rocks. PLoS One 10, e0126941.

545

https://doi.org/10.1371/journal.pone.0126941

546

Xu, S., Payne, M.A., 2009. Modeling elastic properties in carbonate rocks. Lead. Edge 28, 66–74.

547

Yao, Y., Liu, D., Che, Y., Tang, D., Tang, S., Huang, W., 2009. Non-destructive characterization of

548

coal samples from China using microfocus X-ray computed tomography. Int. J. Coal Geol.

549

80, 113–123. https://doi.org/10.1016/J.COAL.2009.08.001

550 551 552 553

Appendix A‒ A‒ Geometric factors 
∗  and 
∗  For dry ellipsoidal pore inclusions, the geometric coefficients are given by the following (Berryman, 1980; Mavko et al., 2009): 
∗  = L MNN , K

(A-1)

26


∗  = O  MNN − L MNN , K

554 555

K

(A-2)

where the tensor MNGB relates the uniform strain field to the strain field within the ellipsoidal

556

inclusion . Berryman (1980) gave the scalar formulation required for calculating 
∗  and 
∗ 

557

as follows: MNN =

LP? P

,

(A-3)

MNN = MNN + K L

Q

P

+

K

PR

+

,

PR PS TPU PV TPW PX P PR

(A-4)

YK = 1 + Z [
\ + ]  − ^  \ + ] − `, L

(A-5)

YK = 1 + Z [
\ + ]  − ^  \ + ] − `, L

(A-6)

L

L

Q

Q

L Q

O Q

L

O

Q

Q

_ _

YQ = 1 + Z [ 1 +
\ + ]  − ^
3\ + 5] ` + a
3 − 4^  + Z
Z + 3a 
3 − L

K

Q

K

Q

Q

4^ \ + ] − ^
\ − ] + 2] Q ,

(A-7)

YL = 1 + Z [1 − \ + ] + ^
\ + ]`,

(A-8)

Y_ = 1 + _ Z\ + 3] − ^
\ − ] ,

(A-9)

L Q

K

YO = Z [−\ + ^ \ + ] − L` + a]
3 − 4^ , _

(A-10)

Yc = 1 + Z1 + \ − ^
\ + ] + a
1 − ]
3 − 4^,

(A-11)

Yd = 2 + Z3\ + 9] − ^
3\ + 5]  + a]
3 − 4^ ,

(A-12)

Yf = Z [1 − 2^ + \
^ − 1 + ]
5^ − 3` + a
1 − ]
3 − 4^,

(A-13)

K _

K Q

K Q

Yg = Z
^ − 1\ − ^]  + a]
3 − 4^ , Z =

"6

"h

(A-14)

− 1,

(A-15)

27

a = L  − " , K 6

^ = ] = \ =

h

K8Qih

Q
K8ih 

l



"6

(A-16)

,

(A-17)

h


 8K

 m
mQ − 1 − cosh8K m ,
m > 1 ?

k   [ cos 8K m − m
1 − mQ  ` ,
m < 1 j
K8 

K8

?

,

(A-18)


3] − 2,

(A-19)

558 559 560 561 562

Table 1. 1. Chemical analysis of the travertine sample. Loss of Ignition (LOI) is the percentage of

563

volatile substances escaping by strongly heating. Since the major concentration goes back to CaO

564

and other oxides are less than 1 %, it can be concluded that the sample is completely composed

565

of calcite. Composition

SiO2

Al2O3

BaO

Value (%)

0.19

0.05

<0.05 55.42 <0.05 <0.05

Composition

MnO

Na2O

P2O5

SO3

TiO2

<0.05 <0.05 <0.05

0.25

<0.05 43.55

Value (%) 566 567 568 569

28

CaO

Fe2O3

K2O

LOI

MgO 0.54

570 571 572 573 574 575 576 577 578 579 580

Table 2. 2. Laboratory measurements of physical properties of core samples measured in Research

581

Institute of Petroleum of Industry (RIPI) and International Geothermal Center (GZB). Length

Diameter

Sample -3

(10 m)

-3

(10 m)

Bulk

Grain

density

density

Weight

VP (m/s)

Porosity

-3

(10 Kg)

VS (m/s)

(%) 3

3

(Kg/m )

(Kg/m )

RIPI

GZB

RIPI

GZB

S1

100.70

54.33

540.99

2320

2720

9.8

5204

5820

2629

2700

S2

100.97

54.48

546.71

2320

2700

9.8

5385

5835

2751

2770

S3

100.37

54.47

566.56

2420

2720

7.4

5530

6012

2811

2990

582 583 584 585 586 587 29

588 589 590 591 592 593 594 595 596 597 Table 3. 3. Properties of mineral and pore phases (Mavko et al., 2009).

598 Phase

Bulk modulus (GPa)

Shear modulus (GPa)

Density (kg/m3)

Calcite

76.8

32

2710

Air

0

0

0

599 600 601 602 603 604 605 606 30

607 608 609 610 611 612 613 614 615 616 Table 4. Results for velocity computations.

617 Sample

Computed VP

Δ % VP-RIPI

Δ % VP-GZB

Computed VS

Δ % VS-RIPI

Δ % VS-GZB

S1

5503

5.7

6.9

2974

13.1

8.7

S2

5436

0.9

8.2

3036

10.3

8.5

S3

5727

3.5

5.8

3087

9.8

2.7

3.4

7.0

-

11.1

6.7

Average 618 619 620 621 622 623

31

624 625 626 627 628 629 630 631 632 633 634

Figure captions:

635

Fig. 1. Core samples used in this study.

636

Fig. 2. 2. An illustration of heterogeneous travertine sample with macro and micro pores.

637

Fig. 3. 3. Multiscale procedure proposed in this study.

638

Fig. 4. 4. Core samples (S1, S2, and S3) and corresponding CT-images extracting as the largest cube.

639

Pore structures obtained via bimodal thresholding with a rough threshold value.

640

Fig. 5. RVE analysis for three samples. A plateau is observed above 1503 voxels for almost all

641

samples.

642

Fig. 6. Computation of porosity cube by three-phase segmentation of the CT-image.

643

Fig. 7. Four examples of 2D high-resolution microscopic and corresponding segmented images.

644

Fig. 8. Bulk and shear moduli of 2D images (circles) and corresponding 3D models (square).

645

Micro-trends based on modified DEM are also shown by dash lines.

32

646

Fig. 9. (a) Lab and computed velocity versus porosity plot and (b) lab versus computed velocity

647

plot.

33

Highlights • • • •

Velocity of a core sample (with 55 mm diameter) is computed via a multiscale procedure. Cost-effective and widely available imaging tools are used. Microscopic 2D images are used to enact 3D effects of pore type in microscale. Velocity computations are highly accurate compared to laboratory tests.

Declaration of interests The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.