Flow Estimation Only from Image Data, based on Persistent Homology

(147 words) 16 Topological data analysis is an emerging concept of data analysis for characterizing shapes. A 17 state-of-the-art tool in topological data analysis is persistent homology, which is expected to 18 summarize quantified topological and geometric features. Although persistent homology is 19 useful for revealing the topological and geometric information, it is difficult to interpret the 20 parameters of persistent homology themselves and difficult to directly relate the parameters to 21 physical properties. In this study, we focus on connectivity and apertures of flow channels 22 detected from persistent homology analysis. We propose a method to estimate permeability in 23 fracture networks from parameters of persistent homology. Synthetic 3D fracture network 24 patterns and their direct flow simulations are used for the validation. The results suggest that 25 the persistent homology can estimate fluid flow in fracture network based on the image data. 26 This method can easily derive the flow phenomena based on the information of the structure.

Abstract (147 words) 16 Topological data analysis is an emerging concept of data analysis for characterizing shapes. A 17 state-of-the-art tool in topological data analysis is persistent homology, which is expected to 18 summarize quantified topological and geometric features. Although persistent homology is 19 useful for revealing the topological and geometric information, it is difficult to interpret the 20 parameters of persistent homology themselves and difficult to directly relate the parameters to 21 physical properties. In this study, we focus on connectivity and apertures of flow channels 22 detected from persistent homology analysis. We propose a method to estimate permeability in 23 fracture networks from parameters of persistent homology. Synthetic 3D fracture network 24 patterns and their direct flow simulations are used for the validation. The results suggest that 25 the persistent homology can estimate fluid flow in fracture network based on the image data. 26 This method can easily derive the flow phenomena based on the information of the structure. 27 28 Introduction 29 Fluid flow processes are ubiquitous in the world, and most are governed by the 30 geometry and nature of the surrounding structures. In particular, recent miniaturization of 31 artificial devices has led to the need for understanding and controlling flow in finer structures. 32 It is also attracting attention to understand flow behaviors in complex fracture networks in 33 developments of natural resources, as in the case of shale gas and geothermal developments. 34 It has been a long-term scientific challenge to predict flow behavior of porous media 35 from structural properties. Permeability is a key parameter for examining flow phenomena in 36 porous media 1 . Permeability cannot be determined only from structure data, and needs to be 37 obtained from laboratory experiments or numerical fluid flow simulations. In contrast, porosity 38 is a parameter that is often used to characterize the structures. The porosity-permeability 39 correlation has been studied extensively in the literature to estimate permeability using porosity 40 (so-called Kozeny-Carman equation) 2,3 . This Kozeny-Carman equation provides a relationship 41 between structure and flow. The correlation has been studied extensively in the literature to 42 estimate permeability using porosity. However, no matter how much void there are, if they are 43 not connected, water cannot flow. Therefore, the Kozeny-Carman equation does not always 44 work. The correlation has been modified to represent real phenomena by adding parameters 45 such as fractal dimension, and tortuosity 2 . These additional parameters can only be determined 46 by fitting, which is not the best way to go about flow prediction based on structural information. 47 Let us also consider flow in a channel from an inlet to an outlet. Hagen-Poiseuille 48 equation is a physical law that describe a steady laminar flow of a viscous, incompressible, and 49 Newtonian fluid through a circular tube of constant radius, r. This is an exact solution for the 50 flow, can be derived from the (Navier-) Stokes equations, and is another way of expressing the 51 relationship between structure and flow.
Since the flow rate is proportional to the cube of the fracture aperture, this relationship between 63 flow and aperture is well-known as the "cubic law" 4-7 . 64 Eqs.
(1) and (2) are only ways to obtain simplified analytic solution to describe the 65 relationship between the flow and structures. This is another way to predict permeability from 66 structural properties 8 . In natural rocks, it is not always a single fracture, but multiple fractures 67 that form a network. Thus, it is necessary to understand not just an individual fracture, but how 68 channels are connected from an inlet to an outlet in whole networks. There have been many 69 studies focusing on networks, but most of the parameters used to describe the structure are 70 probabilistic variables that capture individual fractures, and no suitable parameters have been 71 found yet to evaluate the flow of the entire networks 9 . 72 Topology, a branch of modern mathematics, is good at roughly investigating the 73 connectivity of shapes. Topology focuses on the properties (called topological properties or 74 topological invariants) that are preserved when some form (shape or space) is continuously 75 deformed (stretched or bent, but not cut or pasted). Topology can extract global features that 76 are difficult to capture with machine learning and convolutional neural networks, so it is 77 promising as a complementary feature to extract image information that cannot be detected 78 with other methods. It can be applied to volumetric data as well, so it can pick up information 79 that has been missed in one-way slice-by-slice analysis common to many forms of data 80 processing. 81 Several studies used topological invariants to describe pore-scale structures in porous 82 materials and fracture networks 10,11 . The Minkowski functionals can be interpreted as area, 83 perimeter, or the Euler characteristic, which is a topological constant and were used to link to 84 hydraulic properties 12,13 . Scholz et al. (2019) 16 reviewed the theoretical basis of the Minkowski functionals and its application to 88 characterize porous media. Counting the number of holes using topological invariants like they 89 did is a clue to the shape of the object, and the "essential information" can be extracted well. 90 On the other hand, topology too narrowly focuses on the essential information, it also discards 91 a lot of information, such as size of pore space. The size information, such as radii of tube or 92 apertures of fractures in Eqs.
(1) and (2), must be detected to determine permeability derived 93 analytically. Therefore, previous studies had to add the size information in other ways. 94 Homology is a standard technique for identifying a topological space. In particular, the 95 concept of homology has traditionally played a role in feature extraction focusing on the 96 existence of "holes". Here, the "hole" structure can be regarded as a connected flow channels 97 from an inlet to an outlet. It is expected that topology can be used to detect such connected 98 flow channels. 99 By tracking the sequence of topological spaces, namely, by recording how long 100 homological features persist, we can add information about the size and length of the holes. 101 This can give us a quantitative indication of the size of the holes and the amount of space 102 available, which is called persistent homology. Persistent homology is one of the most 103 important tools in topological data analysis and is expected to compute geometric and 104 topological features of various shapes with ease of computation 17-20 . Thus, this has been 105 applied in several research fields 21-25 , and is also beginning to be used in the analysis of porous 106 materials 26-31 . 107 At this point, in contrast to topological invariants, persistent homology can provide a 108 lot of information that we might need, but it is difficult to interpret the parameters of persistent 109 homology themselves 27-30 . Ushizima et al. (2012) 26 estimates permeability of porous rocks by 110 using Reeb graphs to represent the pore networks. They use persistent homology to distinguish 111 between significant and "noisy" pore spaces, and to supplement the Reeb graphs. Their paper 112 did not go into quantitative evaluation, but focused on qualitative evaluation and visualization. 113 As mentioned before, the "hole" structure that is characterized by topology, can be regarded as 114 connected flow channels from an inlet to an outlet. The aim of our study is to detect the flow 115 channels by persistent homology. Suzuki et al. (2020) 31 proposes a method to detect flow 116 channels in 2D images from persistent homology through image processing. By using their 117 image processing procedure, persistent homology is expected to detect such connected flow 118 channels in complex fracture networks and would also provide their size information such as 119 apertures to predict the permeability. 120 In this study, we applied persistent homology to estimate permeability in fracture 121 networks based on image data. Persistent homology was used to detect the number of flow 122 channels and their apertures in the networks. Synthetic fracture networks were generated, and 123 direct flow simulation was conducted. Permeability derived from persistent homology and 124 simulation results were compared. We applied the proposed method to several published image 125 data and discussed the applicability of permeability estimation based on persistent homology. 126 127 Results

128
Detection of flow channels from persistent homology analysis 129 An example of a fractured rock model with a flow channel connecting an inlet to an 130 outlet is shown in Figure 1a. The yellow area is a solid skeleton, while the white area is 131 fractures forming void spaces. The connecting fractures from the top (inlet) to the bottom 132 (outlet) can be a flow channel. In persistent homology analysis, such structure is recognized as 133 "hole" and quantified as a 1-dimensional hole. Additionally, a discrete island (i.e., connected 134 component) is quantified as a 0-dimensional hole, and a ball (i.e., enclosed solid voids) is 135 quantified as a 2-dimensional hole. The numbers of k-dimensional holes (the dimension of the 136 kth homology vector space) are known to the kth Betti number (β0, β1, and β2). This study 137 focuses on "hole" structures penetrating from an inlet to an outlet, which can be flow channels 138 hence we only analyze 1-dimensional holes in this study. 139 One aspect of persistent homology analysis is that independent fractures are recognized 140 as 2-dimensional holes, as shown in Figure 1b. These independent fractures would not 141 contribute to the fluid flow. Therefore, we can distinguish fractures that act as flow channels 142 and independent fractures by 1-dimensional holes and 2-dimensional holes. arrow) and an isolated pore (mentioned as green). c Schematic of filtration process. 148 Image at t = -1 is the original image. White grids express void spaces. Blue grids are 149 the grids that were removed during the thinning process. Red grids are the grids that 150 were added during the thickening process. represented as sphere cloud data using a pore-network extraction method, then analyzed 158 by calculating the Vietoris-Rips complex topology of the input sphere cloud data. We used an 159 open software HomCloud (https://homcloud.dev/) to analyze binarized 3D images, which can 160 obtain the information of persistent homology by calculating the Euclidian distance of 2D or 161 3D black and white images 32 . 162 Figure 1c shows an example of data process in our persistent homology analysis, 163 called filtration 17 . In filtration, the solid skeletons (yellow parts) are made thinner or thicker, 164 voxel-by-voxel. The original image is set to -1. The process of thinning yellow voxels adjacent 165 to white voxels is regarded as -1, while the process of thickening yellow voxels adjacent to 166 white voxels is regarded as +1. When we reduce the time, the space eventually becomes empty. 167 The nested sequence of the topological spaces from the empty space to the filled space is 168 recorded. The times when the hole appears or disappears are called "birth time" or "death time", 169 expressed as "b". or "d", respectively. 170 In filtration (Figure 1c), a point of the flow channel is closed at t = 1. This closed point 171 is the narrowest aperture in the flow channel. Taking advantage of this, the length of narrowest 172 aperture can be obtained as death time d multiplied by two and its resolution δ. (narrowest 173 aperture = 2 × d × δ = 2 × 1 × 5 = 10). It has been known that the narrowest width in flow 174 channels, which is called critical pore radius, correlates with permeability better than other pore 175 radii 35-37 . Detecting narrowest aperture by persistent homology can therefore be useful to 176 estimate the permeability. 177 The set of pairs (bi, di) for k-dimensional holes is called k th persistence diagram, PDk. 178 If Here is something to keep in mind. A 1-dimensional hole detected by persistent 185 homology is a flow channel penetrating from an inlet to an outlet. At the same time, a ring-186 shaped, internal void-structures is also detected as a 1-dimensional hole. Figure 2a shows an 187 example of a ring-shaped internal void structure. This structure does not connect to the outside 188 (i.e., no flow channel). However, during filtration, the internal void space is closed at t = 1 (d 189 = 1). If images include such hole structures, it would overestimate the number of flow channels. 190 Now, let us prepare an inverted image that the yellow and white are reversed as shown in Figure  191 2b. In filtration, a ring appears at t -1, and the ring width is detected by the value of b (b = -1). 192 Therefore ,  Another example is shown in Figure 2c. This is a ring-shaped internal void structure 198 with two channels that are connected to the outside. In this case, there are two 1-dimensional 199 holes (β1 = 2) with d = 1 and d = 2. Figure  void-structure that is not connected to the outside, and b its inverted image. c ring-209 shaped, internal void-structure with two channels that is connected to the outside and 210 forms a flow channel. d its inverted image. The left column shows 3D view of images.

211
The center column describes processes of filtration. The right column lists Betti 212 numbers ! . 213 214

Synthetic fracture network 215
Synthetic fracture networks were generated by using OpenSCAD 216 (https://www.openscad.org/). We distributed multiple penny-shaped fractures by controlling 217 the apertures, radii, numbers, and orientations of fractures to generate a fracture network 38 . By 218 hollowing out the generated fracture network from a rectangular block, a fractured model 219 where the void spaces were composed of the fracture network was created. This study 220 characterizes one-dimensional flow. The top surface was an inlet, and the bottom surface was 221 an outlet. The fractures were connected from the top to the bottom surfaces. The side 222 boundaries were closed, and water did not flow out from the side. 223 The fractured model is shown in Figure 3. Figure 3a and 3b are the outside and the 224 inside of the model. The fracture orientation was either orthogonal or random. The orthogonal 225 models distributed perpendicular or horizontal fractures to the flow direction (Figure 3c), while 226 the random models distributed fractures by random numbers (Figure 3d). We prepared nine 227 orthogonal models and seven random models. The model parameters for each model are listed 228 in Table 1.  resolution δ of 0.1 mm) were binarized and analyzed by persistent homology using 239 HomCloud 32 . The image size was 360 × 360 × 500 voxels. 240 The Derivation of permeability 256 We use Eq.
(2) to derive permeability, which can be originally calculated by comparing 257 the Stokes equation with the Darcy's law. If Assuming that a fracture is a smooth and parallel 258 plate with the aperture of h and that there is a uniform pressure gradient in one direction within 259 the plane of the fracture, the total volumetric flowrate in the fracture can be written as 260 261 where w is the width of the fracture, perpendicular to the flow direction. h is the aperture, and 264 μ is the water viscosity, dP/dx is the pressure gradient. Darcy's law describes one-dimensional 265 fluid flow through porous media as 266 267 where A is the cross-sectional area. Comparison of Eqs. (3) and (4)  There is an unknown parameter wi in Eq. (6). The 3D voxel data can be regarded as a series of 283 2D cross-sectional images. The 2D cross-sectional image data provides total area of pore space, 284 Ap in each layer. If we introduce effective depth . that is the same for all flow channels, . Estimation of permeability from persistent homology analysis 292 Before applying complex fracture networks, we validated Eq. (7) and our simulation 293 with simple fracture models. Simple models with one or two fractures penetrating from an 294 inlet to an outlet were used (see Figure 5a). Apertures and number of fractures in each model 295 are listed in Table 2. Direct flow simulation with the same fracture network was conducted in 296 OpenFOAM (https://www.openfoam.com/). We could obtain volumetric flow rate and 297 pressure gradient between the inlet and the outlet to calculate equivalent permeability based on 298 Darcy's law. Comparison of permeability between flow simulation and persistent homology 299 analysis is shown in Figure 5b, and listed in Table 2. The persistent homology estimation is in 300 very good agreement with the simulation results. The results show that for such a simple system, 301 persistent homology can estimate the permeability well using Eq. (7). 302 303 Estimation from PH analysis S1 1.0 1 10 1 1.80 × 10 -9 1.80 × 10 -9 S2 1.0 and 1.0 2 10 and 10 2 3.66 × 10 -9 3.47 × 10 -9 S3 0.5 1 5 1 2.39 × 10 -10 2.17 × 10 -10 S4 0.5 and 0.5 2 5 and 5 2 4.79 × 10 -10 4.33 × 10 -10 S5 0.1 1 1 1 1.95 × 10 -12 1.73 × 10 -12 S6 0.1 and 0.1 2 1 and 1 2 3.89 × 10 -12 3.47 × 10 -12 S7 0.5 and 1.0 2 5 and 10 2 2.00 × 10 -9 1.95 × 10 -9 between persistent homology (PH) analysis and direct simulation. The calculated 309 permeability is listed in Table 2. 310 Next, we applied Eq. (7) to complex fracture networks listed in Table 1. Comparison of  311 permeability between flow simulation and persistent homology analysis is shown in Figure 6. 312 The estimation is in reasonable agreement with the simulation results. 313 314

Figure 6 | Estimation of permeability by persistent homology (PH) analysis for 315
orthogonal fracture networks (blue) and random fracture networks (orange). The 316 calculated permeability is listed in Table 1.
There is a limitation of Eq. (7). Eq. (7) is based on a parallel plate model, so the flow 319 is assumed to be straight. If there is tortuosity in a flow channel, the flow length will be longer, 320 and the estimated permeability may be larger than the true value. Figure 7 shows streamlines 321 in model O8 colored as green. We can see that the streamlines are winding and flowing. Keep 322 in mind the fact that tortuosity was not taken into account in Eq. (7). 323   1.54 × 10 -9 GL-D1 1.92 × 10 -10 2.74 × 10 -10 GL-D2 1.70 × 10 -10 4.96 × 10 -10 GL-D3 1.47 × 10 -10 7.98 × 10 -10 GL-D4 1.44 × 10 -10 9.71 × 10 -10 GS-D1 9.77 × 10 -10 5.19 × 10 -9 GS-D2 9.75 × 10 -10 3.34 × 10 -9 GS-D3 9.60 × 10 -10 2.43 × 10 -9 GS-D4 9.18 × 10 -10 2.67 × 10 -9 P-D1 3.25 × 10 -10 1.75 × 10 -9 P-D4 4.01 × 10 -10 1.53 × 10 -9 Berea 2900 x 2320 to represent the pore networks. They use persistent homology to distinguish between 394 significant and noisy pore spaces, and to supplement the Reeb graphs. In fact, the Reeb graph 395 and persistent homology were used independently and separately. We think that using Reeb 396 graphs is a good direction to go to the next step. 397 We have succeeded in modeling physical phenomena from image data based on the 398 topological data analysis. The method could be applied also to a wide range of porous media 399 including artificial devices. It is expected to be applicable not only to estimate flow properties, 400 but also to characterize different transport phenomena, such as mass transfer, electrical and 401 magnetic flows. 402 403 Method 404 Persistent homology analysis 405 STL files of synthetic fracture networks were generated by using OpenSCAD, and the 406 STL files were converted to the cross-sectional images in 36 mm × 36 mm × 50 mm with a 407 voxel resolution of 0.1 mm in Autodesk Netfabb. To eliminate some unexpected noises, all the 408 images were blurred in XnConvert. The png files of the image data were analyzed in 409 HomCloud 63 . When there were small differences between birth and death of PH1, the hole 410 structure may appear during the image analysis. Thus, we neglected the result with di -bi < 2 411 were eliminated. boundaries. b discretized model.