Investigating rough single-fracture permeabilities with persistent homology

The permeability of rock fractures is a crucial parameter for flow processes in the subsurface. In the last few decades, different methods were developed to investigate on permeability in fractures, such as flow-through experiments, numerical flow simulations, or empirical equations. In recent years, the topological method of persistent homology was also used to estimate the permeability of fracture networks and porous rocks but not for rough single fractures yet. Hence, we apply persistent homology analysis on a decimetre-scale, rough sandstone bedding joint. To investigate the influence of roughness, three different data sets are created to perform the analysis: (1) 200 µm, (2) 100 µm, and (3) 50 µm resolutions. All estimated permeabilities were then compared to values derived by experimental air permeameter measurements and numerical flow simulation. The results reveal that persistent homology analysis is able to estimate the permeability of a single fracture, even if it tends to slightly overestimate permeabilities compared to conventional methods. Previous studies using porous media showed the same overestimation trend. Furthermore, the expenditure of time for persistent homology analysis, as well as air permeameter measurements and numerical flow simulation, was compared, which showed that persistent homology analysis can be also an acceptable alternative method.

1 Introduction

The permeability of rocks is a crucial parameter for fluid flow processes in the subsurface. While the prevailing flow processes in porous media are well understood, a different picture emerges when the flow is dominated by fractures (Suzuki et al., 2019). Although flow in geological settings controlled by fractures is occurring in both shallow aquifers (e.g. for drinking water supply) and deep reservoirs (e.g. geothermal energy production and oil and gas abstraction), fractures are often simplified as two parallel plates. In addition, due to complexity reasons, a single parameter, the hydraulic aperture, is often used to represent the permeability of a single fracture or even entire discrete fracture networks (Min et al., 2004; Blum et al., 2009; Müller et al., 2010). However, due to roughness of fracture surfaces, a single value is not sufficient to capture the flow channelling and critical flow paths (Tsang, 1992; Tsang and Neretnieks, 1998). Hence, investigations of the fracture roughness are essential, although this is more expensive in terms of the costs and time (Tatone and Grasselli, 2012, 2013).

Nowadays, various methods are available to study how fluid flow is influenced by fracture roughness through single fractures. These can be divided into four major categories: (1) empirical, (2) experimental, (3) numerical, and (4) geometric methods.

Empirical methods are simple and also fast, cheap, and often sufficient to provide a first overview over the flow behaviour of a fracture. There are different empirical models derived from flow experiments, numerical simulations, or statistical models of different fracture types (Louis, 1972; Barton and Quadros, 1997; Xiong et al., 2011; Kling et al., 2017; Suzuki et al., 2017). Often, solely mechanical apertures and relative roughness, depending on the standard deviation, are required to not only calculate hydraulic apertures but also the fractal dimension or the Peklenik number, which is defined as the ratio of the correlation length in the x and y directions (Patir and Cheng, 1978; Brown, 1987).

Experimental methods provide more detailed results than empirical methods which use more simplified relations for fast practical application. The standard methodology is flow-through experiments at laboratory scale to observe flow patterns in single fractures (Brown et al., 1998; Watanabe et al., 2008; Ferer et al., 2011). In recent years, hydromechanical coupling (Vogler et al., 2018; Wang et al., 2021) or reactive transport have additionally been performed for the exclusive investigation of flow patterns in flow-through experiments (Durham et al., 2001; Huerta et al., 2013). In addition to typical flow tests at laboratory scale, it is also possible to perform flow experiments on larger scales in the laboratory or in the field (Novakowski and Lapcevic, 1994; Thörn and Fransson, 2015; Weede and Hötzl, 2005). Flow-through experiments allow us to investigate direct fluid flow through fractures. The flow distribution and preferred flow paths can be predicted by replicating the fracture geometry in transparent materials. Beside flow-through experiments, air permeameters can be used to determine the permeability of fractures (Cheng et al., 2020; Hale et al., 2020; Hale and Blum, 2022). Air permeameters allow us to measure the permeability of fractures directly on an outcrop or drilling core (Brown and Smith, 2013). In addition, it is also possible to obtain a zonal observation of the permeability, since several measurements have to be conducted along a fracture outcrop due to the measurement method (Hale et al., 2020).

Another way to investigate flow in fractures is the representation of the fracture in numerical models. For this purpose, the geometry of a fracture is either projected onto a two-dimensional surface (Pyrak-Nolte and Morris, 2000; Javanmard et al., 2021) or represented in three dimensions (Javadi et al., 2010; Xiong et al., 2011; Wang et al., 2016; Chen et al., 2021). The latter increases the computational effort considerably and also allows a more accurate investigation of flow processes. Another major advantage is that it is possible to simulate various scenarios under different conditions, such as high confining pressures or high flow rates, which would exceed technically possible conditions in laboratory experiments. In addition, numerical models are able to consider flow effects, region by region, in the fracture and therefore characterise the main flow paths (Marchand et al., 2020; Javanmard et al., 2021). However, the geometry of the fracture and the prevailing boundary conditions for simulations have to be precisely known in order to obtain meaningful results (Barton et al., 1985; Tatone and Grasselli, 2012).

If the focus is set on the investigation of the flow behaviour and permeability distribution within the fracture, then geometric methods can also be used. The Kozeny–Carman equation is well known as a representative method for estimating flow properties from pore structures (Kozeny, 1927; Carman, 1937), and attempts have long been made to estimate permeability by extracting porosity from images without experiments or numerical simulations (Costa, 2006; Torskaya et al., 2014; Oliveira et al., 2020). More recently, attempts have been made to use machine learning or deep learning on images (Sudakov et al., 2019; Anderson et al., 2020; Araya-Polo et al., 2020; Hong and Liu, 2020; Alqahtani et al., 2021; Da Wang et al., 2021). Of course, deep learning is a powerful method, but adequate use of machine learning requires deep technical understanding, rigorous testing, and sufficient amounts of training data.

Topological data analysis (TDA) is another way to extract the information of shapes and structures from big data (Carlsson, 2009; Thiele et al., 2016). TDA is an analysis method that focuses on the structure of data within the field of algebraic topology, demonstrating particular strengths in handling data types such as images, complex structures, and networks. TDA can capture the structure of the data in a rough sense and characterise the features of connectivity and holes, therefore ignoring noises in data and extracting important information, such as the independent of coordinate system and number of dimensions. Persistent homology (PH), one of the most used TDA methods, can capture changes and continuity of topological features by tracking algebraic descriptions called homology. This method was developed in the early 2000s (Edelsbrunner et al., 2000; Zomorodian and Carlsson, 2005) and is applied in various research fields such as materials science (Hiraoka et al., 2016), computer science (Choudhury et al., 2012), and biology (Chan et al., 2013). In geosciences, this approach has only been applied in the past decade to characterise porous rocks and to determine their permeability (Delgado-Friedrichs et al., 2014; Robins et al., 2016; Bizhani and Haeri Ardakani, 2021). Furthermore, the determination of hydroelastic properties of porous media is possible with this method (Jiang et al., 2018). In the field of fractured rocks, persistent homology was recently also applied to study small-scale fracture networks (Suzuki et al., 2020, 2021). Based on these studies, the general application of PH for permeability estimation of fracture networks could be demonstrated. In these small-scale (millimetre to centimetre scale) studies, the effect of fracture roughness was not particularly investigated (Suzuki et al., 2020, 2021). Further research is therefore needed to investigate larger-scale fractured rocks, in which surface roughness has a significant effect on flow behaviour.

The objective of this study is therefore the application of the persistent homology analysis on a natural, mesoscale (decimetre scale) single fracture to estimate the permeability. The focus is on the anisotropy of permeability in different flow directions, as well as the influence of resolution on permeability. Thus, three data sets of the same fracture are prepared, which have different resolutions (50, 100, and 200 µm). Finally, these results are compared with results from experimental air permeameter measurements, as well as numerical flow simulations.

2 Methods

2.1 Fracture sample

The fracture sample is a natural bedding joint in a sandstone block taken from a quarry in Bebertal, Germany (Fig. 1; Heidsiek et al., 2020; Hale and Blum, 2022). The sandstone is Flechtinger sandstone, an oil and gas reservoir rock in the North German Basin. The block contains one bedding joint, with a length of 120 mm in the x direction and 450 mm in the y direction (Fig. 1). Previous studies have already characterised the relevant hydromechanical properties of Flechtinger sandstone such as porosity (9 %–11 %), matrix and fracture permeability, Young's modulus and bulk modulus, thermal dependencies of stress and strain behaviour, and the mineralogical composition (Frank et al., 2020; Cheng et al., 2020; Fischer et al., 2012; Hale and Blum, 2022; Heidsiek et al., 2020; Hale et al., 2020; Hassanzadegan et al., 2012; Blöcher et al., 2019). Of particular interest for this study is the low-matrix permeability of 0.1–1 mD (Cheng et al., 2020; Hassanzadegan et al., 2012). Furthermore, the findings of Gutjahr et al. (2022), Hale et al. (2020), and Hale and Blum (2022) are crucial, since they performed investigations on fracture permeability on exactly the same fracture. Gutjahr et al. (2022) investigated on the roughness of the fracture and calculated the Hurst exponent for different angles. The medians of all Hurst exponents in the x direction and in the y direction are 0.48 and 0.42, respectively. Hale et al. (2020) and Hale and Blum (2022) determined the average fracture permeability to be $\mathrm{5.6}\times {\mathrm{10}}^{-\mathrm{10}}$ m². In addition, they found that the centre of the fracture is less permeable than the left and right side of the fracture (according to the front view of the sandstone block shown in Fig. 1a). On the right-hand side of the fracture, this can be explained by a barite vein intersecting the fracture, which was formed before the fracture opening. In the closer vicinity of the vein, the mechanical aperture is increased compared to central parts of the fracture. Comparing fracture and matrix permeabilities shows that the fracture permeability exceeds the matrix permeability by more than 8 orders of magnitude. Thus, the matrix permeability is considered negligible in this study.

Figure 1 (a) Photo of the studied sandstone block showing also the investigated bedding joint. The front surface of the real block corresponds to the y–z plane in the digital model on the right-hand side. (b) 3D model of the bedding joint and the surrounding sandstone block.

2.2 Data preparation and matching

For data preparation, a self-developed Python code called FracMatchPy was used. This code is able to (1) match two separated and spatially uncorrelated fracture surfaces and (2) create binary cross sections of the fracture as input for the PH analysis.

In the matching process, meshed laser scans of the bedding joint surfaces were used as input data for the Python code (Fig. 2). The surfaces were scanned by a combined system of the high-resolution laser scanner Nikon ModelMaker MMDx 100 and the articulated arm MCA II on which the scanner was mounted (Nikon Metrology NV, 2010, 2018). The scanner provides a resolution of 100 µm and a single-point accuracy of 10 µm (Nikon Metrology NV, 2018). The scanned point cloud was then meshed using MeshLab (Cignoni et al., 2008). Since the meshes were not spatially related, the two surfaces were roughly matched by hand to shorten the runtime of the Python code for data preparation. The exact matching was then performed stepwise by several rotation and translation steps within specified limits using the Python code. The translation limit was 3 mm total distance in each direction so that one surface was displaced in 100 µm steps, starting from the geometric centre of the surface. In addition, there was a rotation range of −0.3 to 0.3° in which the rotation was performed in 0.01° steps. The best fit was then determined by applying a minimisation function of the average mechanical aperture using the arithmetic mean between the fracture surfaces.

Figure 2 Workflow of this study with three distinct steps, which is applied on the laser-scanned fracture surfaces.

2.3 Binary image generation and permeability estimation using persistent homology

We apply the permeability estimation method proposed by Suzuki et al. (2021), which uses persistent homology (PH) to extract information about the flow channels from image data. The first step is to convert the fracture information prepared in Sect. 2.2 into 3D binarised image data sets.

Three-dimensional image construction of the bedding joint and the surrounding sandstone block was generated. The image construction is a series of binary cross-sectional images of the x–z planes along the y axis. The size of the image area is 115 mm × 8.4 mm × 392.2 mm. Three data sets with different resolutions were generated to investigate the effect of varying resolution on permeability estimation. The first data set was created with a low resolution of 200 µm in all spatial directions. The data set contains $\mathrm{578}\times \mathrm{1962}\times \mathrm{42}$ voxels. The second data set was generated with a medium resolution of 100 µm in each spatial direction. The data set contains $\mathrm{1154}\times \mathrm{3922}\times \mathrm{84}$ voxels. In the third data set, the resolution was again reduced by half to 50 µm in all spatial directions. The data set contains $\mathrm{2308}\times \mathrm{7844}\times \mathrm{169}$ voxels. The fracture was considered to be fully permeable, with no low permeable filling or sealing. Since the matrix permeability is 8 orders of magnitude lower than the fracture permeability, the matrix was considered to be fully impermeable (Cheng et al., 2020; Hassanzadegan et al., 2012). Thus, each image contains the permeable fracture in white colours (binary value = 1) and the impermeable matrix in black colours (binary value = 0).

PH analysis was then performed using HomCloud (Obayashi et al., 2022). Since HomCloud can only handle a maximum data volume of 1021³ voxels in a single-analysis run, the data sets were divided into several subpackages which were processed separately (Fig. 2). The data set with low resolution (200 µm/pixel) has $\mathrm{578}\times \mathrm{1962}\times \mathrm{42}$ voxels and was split into two divisions in the y direction. Thus, two subpackages were created for the data set with low resolution, each of which contained $\mathrm{578}\times \mathrm{981}\times \mathrm{42}$ voxels. Dividing the x direction and z direction for the data was not necessary because the length of one side of the image was <1021 pixels. The data set with middle resolution (100 µm pixel) has $\mathrm{1154}\times \mathrm{3922}\times \mathrm{84}$ voxels and was split into two divisions in the x direction and four divisions in the y direction. The data set was divided into eight subpackages, with an average of $\mathrm{577}\times \mathrm{981}\times \mathrm{84}$ voxels each. The data set with high resolution (50 µm/pixel) has $\mathrm{1154}\times \mathrm{3922}\times \mathrm{84}$ voxels and was split into three divisions in the x direction and eight divisions in the y direction. The data set was divided into 24 subpackages, with an average of $\mathrm{769}\times \mathrm{981}\times \mathrm{169}$ voxels each.

PH analysis was then performed for each subpackage. HomCloud enables the extraction o topological features; a 3D object can consist of three different topological features. First, 0-dimensional topological features characterise connected components such as impermeable, solid matrices without voids. Structures such as fractures or connected pores are recognised as 1-dimensional (1D) topological features. Finally, 2-dimensional (2D) topological features are represented by enclosed pores not connecting to other flow channels. Since only fractures connected between an inlet and an outlet of the domain can serve as potential flow channels, this study focuses on the analysis of 1D topological features.

In HomCloud, a process called filtration is performed to detect 1D topological features from image data. For simplicity, an example of a two-dimensional image is shown in Fig. 3. In the binary image, the black areas are the rock matrix (binary value = 0), while the white areas are the void space (binary value = 1). During the filtration, the black pixels are thinned or thickened by 1 pixel from the boundary between the white and black pixels. The process is considered to be time variation, and the initial image is assumed to be time 0 (t=0). The time change in the negative direction is to make the black pixels thin, and the time change in the positive direction is to make the black pixels thick. The images are stored at each time. Continuing this operation, if the pixels continue to be thinned (time change in the negative direction), then the image becomes all white. If time is advanced in the positive direction from the all-white state, then a channel (i.e. 1D topological feature) appears at certain times, and if time is continuously advanced in the positive direction, then the channel disappears further. Note that the channel (1D topological feature) here is the connected shape from left to right, surrounded by black pixels. These times are called birth (b) and death (d), respectively. This method is named “persistent homology” because it attempts to see how persistent topological features are. While other topological data analysis extracts only topological features, the advantage of PH is that by utilising birth–death information, one can obtain not only topological features but also the geometric information on length.

In HomCloud, the output is a persistence diagram representing a frequency distribution of the number of birth and death pairs, as shown in step II in Fig. 2. From the definitions of birth and death, the presence of a flow channel (1D topological feature) in an image means that birth time is negative and death time is positive. Thus, the number of such pairs ($b<\mathrm{0}<d$) can be considered to be the number of flow channels in the image. In addition, in the case of the process of the thickening of the black area, as shown in Fig. 3, the fracture is closed from both sides (see the image labelled “t=d”). Therefore, the doubling of death and multiplying by the resolution of the image can be taken as the smallest aperture of the channel. Thus, from HomCloud, the frequency distribution of the number of channels present in the image and their minimum aperture widths can be obtained. It is important to note that the 1D topological features evaluated in PH include the aforementioned left-to-right connected shapes surrounded by black pixels, as well as the void ring structures that are not connected to the outside. How to remove such structures is described in Suzuki et al. (2021). It should also be noted that if the channels are connected like a ladder, then PH may detect a large number of channels.

Suzuki et al. (2021) assumed that the channels have parallel-plate geometries and that they are parallel to each other. The permeability is estimated based on the power law as follows:

$$\begin{array}{}\text{(1)}& K=\sum _{i=\mathrm{1}}^N{\displaystyle \frac{w_i{h}_i^{\mathrm{3}}}{\mathrm{12}A}}.\end{array}$$

A is the surface area of the cross section of the medium, and N is the number of flow channels. w_i is the depth of flow channel i, and h_i is the aperture of the flow channel i. As mentioned above, the number of flow channels N was estimated from the number of birth–death pairs, and the aperture h_i was estimated as 2d_iδ in PH analysis, where d_i is the death of the flow channel i, and δ is the resolution of the images. The average of the depth of flow channel $\overline{w}$ is determined by the cross-sectional area of the image (Suzuki et al., 2021). Thus, the above equation can be converted to the parameters from PH and image analysis.

Figure 3 Schematic illustration of the filtration process during persistent homology analysis. The permeable area of the binary image (fracture) is thickened (cyan) and thinned (pink).

2.4 Experimental measurement of the fracture permeability

An air permeameter was used to experimentally measure the fracture permeability. The transient air permeameter TinyPerm 3 (manufactured by New England Research, Inc.) is a portable device, which is able to directly measure matrix permeabilities or hydraulic fracture apertures on outcrops or cores (Filomena et al., 2014; Cheng et al., 2020; Hale and Blum, 2022). For this purpose, the device is filled with air by lifting the piston, and the rubber nozzle of the instrument is pressed against the fracture outlet. The measurement is started by depressing the piston in direction of the sample, creating a vacuum between the fracture and the device. A microcontroller unit in the air permeameter simultaneously records the transient change in air pressure at the fracture outlet and the volume change in the device. Once the total volume of the device is compressed, the permeability is automatically determined from the recorded curve (Brown and Smith, 2013). The range of measurable permeabilities is from 1 mD to 10 D for porous rocks and hydraulic apertures of approximately 10 µm to 2 mm for fractures (New England Research Inc., 2016). The latter corresponds to permeabilities in orders of magnitude from 10⁻¹¹ to 10⁻⁷ m².

This study uses both existing air permeameter measurements from Hale et al. (2020) and measurements conducted in this study. In the study of Hale et al. (2020), experiments along the long edge (y direction) of the fractured block were performed and the permeability in the x direction was measured. In addition, complementary measurements along the x axis with permeabilities measured in the y direction were performed in the frame of this study. In total, the y axis was divided into 21 sections and the x axis into 4 sections. The corner areas were not included in the measurement due to breakouts from the block. Each section was measured 10 times to obtain average values. The hydraulic aperture along the y edge in the x direction was determined to be 81±1 µm (Hale and Blum, 2022) and that along the x edge in the y direction to be 57±1 µm.

2.5 Numerical flow simulation

Apart from experimental air permeameter measurements, hydraulic apertures of the fracture were also determined by numerical flow simulations, using the Multiphysics Object-Oriented Simulation Environment (MOOSE) framework (Permann et al., 2020). Within this framework, the fluid flow through the fracture was simulated with the SaintBernard application (Schädle, 2020), which is based on the inbuilt PorousFlow module (Wilkins et al., 2020, 2021). In this application, a 3D fracture is projected on a 2D surface embedded in a 3D environment. The aperture of the fracture is assigned as a permeability parameter to each cell of the 2D mesh. The fluid flow velocity is then simulated in the lower dimension, and the hydraulic aperture is calculated, considering Darcy flow and the cubic law with the following equation for each cell of the mesh:

$$\begin{array}{}\text{(2)}& a_{\mathrm{h}}=\sqrt[\mathrm{2}]{{\displaystyle \frac{\mathrm{12}v\mathit{\mu }L}{\mathrm{\Delta }p}}}.\end{array}$$

In Eq. (2), a_h is the hydraulic aperture of the fracture, v is the fluid flow velocity, μ is the dynamic viscosity, L is the length of the fracture (in flow direction), and Δp is the hydraulic pressure gradient. Further information about the numerical simulation and SaintBernard can be found in Javanmard et al. (2021).

Similar to the air permeameter measurements, numerical flow simulations were performed with fluid flow in both the x and y direction to determine the permeability in each of these.

3 Results and discussion

3.1 Permeability estimation from persistent homology analysis

The prepared binary data sets were used to calculate the permeability of the fracture. The permeability was calculated in two different flow directions for flow parallel to the x axis and parallel to the y axis (Fig. 4). The calculated permeabilities for the 200 µm resolution data set are $\mathrm{6.4}\times {\mathrm{10}}^{-\mathrm{10}}$ m² in the x direction and $\mathrm{6.2}\times {\mathrm{10}}^{-\mathrm{10}}$ m² in the y direction. Using the cubic law, the calculated hydraulic apertures are 88 and 86 µm, respectively. For the data set with 100 µm resolution, the permeabilities are $\mathrm{4.4}\times {\mathrm{10}}^{-\mathrm{10}}$ m² in the x direction and $\mathrm{4.0}\times {\mathrm{10}}^{-\mathrm{10}}$ m² in the y direction. This corresponds to hydraulic apertures of 73 µm in the x direction and 69 µm in the y direction. For the 50 µm resolution data set, the permeability in the x direction is $\mathrm{7.0}\times {\mathrm{10}}^{-\mathrm{10}}$ m² and $\mathrm{3.2}\times {\mathrm{10}}^{-\mathrm{10}}$ m² in the y direction, which equals hydraulic apertures of 92 and 62 µm, respectively. Comparing the three different data set shows that PH analysis for all data sets results in higher permeability in the x direction than in the y direction ($k_x/k_y>\mathrm{1.0}$).

The $k_x/k_y$ ratio of the data sets with 200 and 100 µm resolution is nearly identical, 1.0 and 1.1, respectively, whereas the ratio of the 50 µm data set shows a higher ratio of 2.2.

Detailed examination of the individual fracture surfaces and the matched fracture shows that the highest mechanical apertures of >1 mm occur mainly along the barite vein that crosses the fracture parallel to the x direction (Fig. 4). From previous studies on the fracture, it appears that this barite vein dominates the flow behaviour and, thus, forms the main flow path along the fracture due to its increased mechanical aperture and lower roughness compared to other regions of the fracture (Hale et al., 2020). Hence, it serves as a preferential flow path in the x direction, whereas it acts more as a barrier or redirection for flow in the y direction. The study of Gutjahr et al. (2022) shows that, in addition to the anisotropy of the permeabilities, a slight anisotropy of the roughness can be observed. Analogous to the permeability, the Hurst exponent in the x direction is higher (H_x=0.48) than in the y direction (H_y=0.42). It is noteworthy that the ratio of the Hurst exponents ($H_x/H_y$) is 1.1 and thus corresponds well with the determined ratios for permeabilities. This is reasonable, as increased roughness tends to result in more distinct flow channels with larger mechanical apertures. Consequently, this leads to increased permeabilities.

Figure 4 Permeability determined by air permeameter measurements, numerical flow simulations, and persistent homology using data sets with resolutions of 200, 100, and 50 µm.

3.2 Comparison to air permeameter measurements and numerical simulation

In addition to the estimation of the fracture permeability by PH, the results were also compared to permeabilities derived from alternative methods. Thus, a comparison was performed using experimental air permeameter measurements and numerical flow simulations to show the validity of the PH analysis. In Fig. 4, the values for the permeability in the x and y directions, as well as the ratio of the two permeability values of air permeameter measurements and numerical simulation, are also shown. All x permeabilities differ by $<\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{10}}$ m² from the experimental or numerical results. However, the permeabilities in the y direction scatter slightly more. Overall, there is a good fit between PH analysis and alternative methods, which is reflected by a root mean squared error (RMSE) of $\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{10}}$ m². Normalisation with the difference between maximum and minimum of all observed permeabilities leads to a normalised root mean squared error (NRMSE) of 0.34. However, there is also the trend that the values agree increasingly better with the comparative values as the resolution of the data set increases. For example, the RMSE of the x and y permeability between air permeameter or numerical simulation and the 200 µm resolved data set is $\mathrm{2.0}\times {\mathrm{10}}^{-\mathrm{10}}$ m². For the 100 µm data set, it reduces to $\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{10}}$ m², and for the 50 µm data set, it is $\mathrm{1.1}\times {\mathrm{10}}^{-\mathrm{10}}$ m². The NRMSEs are 0.54 (200 µm), 0.38 (100 µm), and 0.25 (50 µm), respectively. This trend was also shown in studies using persistent homology in porous media or fracture networks before (Moon et al., 2019; Suzuki et al., 2021). The study of Moon et al. (2019), in which fluid flow through the pore spaces of different digital sandstone and chalk samples was examined, could show particularly that the number of excessively high outlier permeabilities can be prevented with higher resolution. Similar findings are shown in Suzuki et al. (2021), in which permeabilities that are significantly higher than the comparison simulation could also be reduced by improving the resolution. Conclusively, it can be stated that apart from the permeability estimation of the 50 µm data set in the x direction, the results of the PH analysis improve with increasing resolution.

3.3 Evaluation of persistent homology analysis

Since the results of three different methods for permeability determination are in good agreement (NRMSE = 0.34), a classification of the results in the context of other PH analyses was carried out. The study of Suzuki et al. (2021) applied PH on different data sets of porous and fractured rocks of previous studies (Andrew et al., 2014; Muljadi et al., 2016; Mehmani and Tchelepi, 2017). The 16 data sets of porous media and 15 data sets of fracture networks were analysed, each with PH and numerical flow simulation. In Fig. 5, the results of this study, as well as the results by Suzuki et al. (2021), are shown. The method of PH for permeability estimation in both studies provides comparable results to the respective reference method. The results of this study are closely scattered around the 1:1 line and thus match well with the results based on fracture networks (dark grey diamonds). It appears that the results of this study can be estimated even better than most of the data points of previous data sets, especially those generated from porous media (light grey crosses). In the latter, PH tends to overestimate the permeability, which cannot be confirmed for the data in this study. However, this study indicates that the quality of the permeability estimation is not only attributable to the type of cavity (pores, single fractures, or fracture networks). Based on the results, the quality of the permeability estimation by PH is also dependent on the resolution and anisotropy of the respective data set. Nevertheless, a larger number of data sets should be examined for a more precise assessment of the various influences on the quality of the permeability estimation.

Figure 5 Comparison of the estimated permeabilities of this study with the estimated permeabilities for porous media (light grey crosses) and fracture networks (dark grey diamonds) of previous studies by Suzuki et al. (2021), using data of Andrew et al. (2014), Muljadi et al. (2016), and Mehmani and Tchelepi (2017).

In previous sections, it is shown that PH provides comparable results for permeability estimation of rough single fractures compared to other more conventional methods. However, for it to be an alternative, the effort and computational time also has to be considered. In Table 1, an estimation of different working steps is presented for the permeability estimation of single fractures. For PH analysis, the dependency of the time needed on the resolution is also considered. The table is divided into three working steps. The step of preparation and preprocessing contains every working step after the collection of a sample, including imaging and matching of fracture surfaces, generating binary images, setting up a numerical model, or initialising the air permeameter. In the second working step, the effort for the actual measurement, numerical simulation or PH analysis is considered for one sample. In postprocessing, all working steps recalculating the fracture permeability from the measured, simulated, or analysed results is considered. Comparing the expenditure of time for all methods and resolutions, all methods apart from PH analysis for 50 µm data are ranging in the same order of magnitude. An increase to 50 µm resolution demands an enormous increase in time expenditure, mainly due to the extremely high analysis time. Considering the quality of the results, the data set with 100 µm therefore seems to be an adequate alternative to conventional methods, as it can provide high-quality results with similar efforts.

Table 1 Estimation of expenditure of time in hours for air permeameter measurements, numerical flow simulation, and persistent homology for three different resolution steps (200, 100, and 50 µm).

4 Conclusions

This study shows that persistent homology (PH) provides acceptable results for the permeability estimation of a natural bedding plane joint of a sandstone. Compared to other methods such as experimental air permeameter measurements and numerical flow simulation, it tends to slightly overestimate lower permeabilities. The overestimation of permeabilities is also traceable in previous studies using PH analysis on porous media, as well as small-scale fracture networks. For single-fracture application, a reason could be the application of the cubic law, which tends to overestimate the permeability in fracture networks or rough fractures.

In comparison to other methods, PH is a cheap and time-effective method. Once the geometry of a fracture is imaged (e.g. laser scanning, computed tomography, and structure from motion), all the tools to determine permeability are accessible and open-source. In contrast to most experimental methods, no laboratory is required to estimate the permeability. An exception to this is the air permeameter used here, which is even applicable in field experiments due to its portability. The advantages of PH compared to numerical modelling are, first, the lower required computing capacity and computing time. Furthermore, the number of parameters required to successfully perform a simulation is significantly reduced, since only the geometry is sufficient as an input parameter.

Suzuki et al. (2020) showed that small-scale (millimetre to centimetre scale) discrete fracture networks (DFNs) can be precisely studied by PH analysis. Our study demonstrates the applicability of the methods to a mesoscale (decimetre scale), rough bedding joint of Flechtinger sandstone. In order to verify these results, future work should focus on other types of fractures, such as open-mode or shear-mode fractures, as well as different lithologies, such as fractures in granites or clay. Furthermore, the influence of roughness on the flow behaviour and the permeability distribution across the fracture should be investigated. In addition, a more detailed investigation of the permeabilities of different areas of this fracture could be performed on the basis of the high-resolution scans used here. This could also allow a potential scaling effect of the permeabilities to be analysed in more detail. By combining the approach of fracture networks and rough single fractures, a long-term objective could also investigate the influence of fracture surface roughness on fracture networks at larger scales. In particular, the well-functioning subdivision of the total data set into smaller, high-resolution subpackages should allow analysis of larger DFNs without compromising the resolution. Other potential future work could be found in the comparison of such PH analysis with numerical DFN models. In addition, since numerical models for the accurate representation of fracture networks are computationally expensive, it is expected that PH is able to save time and costs.