Quantitative Analysis of Pore Space Structure in Dry and Wet Soil by Integral Geometry Methods

Abstract

We present a methodology for a numerical analysis of three-dimensional tomographic images in this paper. The methodology is based on integral geometry, topology, and morphological analysis methods. It involves calculating cumulative and non-cumulative pore size distributions of Minkowski functionals and Betti numbers. We investigated 13 samples in dry and wet (saturated beyond the field capacity) conditions within different horizons of the Phaeozem albic. For samples of the arable horizon, an increase in the Euler characteristic was observed in the process of wetting. For samples from the A2, AB and B2 horizons, the Euler-Poincare characteristic decreased during wetting. It has been proven that both Betti numbers (number of isolated pores and number of “tunnels”) decrease with swelling of the AB and B2 horizons at a depth of 20–90 cm. For samples from the arable horizon, another dependence was observed: A Betti number of zero increased first but decreased during wetting. Based on the change in topological characteristics, two methods of changing the topology of the void space of the soil were demonstrated. The above-described quantitative changes of proposed parameters of pore space tomographic images prove the possibility and progressiveness of their usage for the pore space transformation estimate.

1. Introduction

Soil porosity is a three-dimensional structure with complicated internal surface geometry. The mechanical and hydrological soil characteristics are mainly determined by the three-dimensional pore space structure and by a form of phase boundary “pore-solid skeleton” [1,2,3]. X-ray computed tomography allows us to obtain information about pore space structure with a high degree of extension without sample destruction. However, different methods for three-dimensional images quantitative analysis are needed to obtain quantitative information about the pore structure.

There are numerous examples of two-dimensional [4] and three-dimensional [5,6] soil image analysis using statistical, morphological, and other methods. All methods have both advantages and disadvantages, while few have been validated. Methods which have their fundamental basis in statistical physics, integral geometry, and topology are examined in the present paper. These methods are based on Minkowski functionals, i.e., four basic geometric and topological characteristics that are associated with three-dimensional objects with the purpose of structure description.

Investigations have been devoted to the application of methods based on integral geometry and morphological analysis. These methods are used for analysis media with complex internal structures at different scales [7] including soils, sedimentary rocks, foams [8], ceramics [9], and composite materials [10]. Early studies focused on artificial media [11,12] and 2-dimensional rock samples cuts [13]. In later works, 3-dimensional X-ray computed tomographic [14,15,16] and FIB-SEM [17] images were examined. The effect of the resolution of 3-dimensional images on pore-size distribution and Minkowski functional values were discussed in [18].

Specifically, in [16], intrinsic geometry quantitative analysis was examined. The samples were managed in several ways based on Minkowski functionals evolution during the process of dilatation either of the solid phase or voids analysis. This method is applicable to objects at all scales—the integral geometry and topological methods can quantitatively describe the geometry of voids in soil aggregates [19] or voids between aggregates or macroscopic pores [20].

We would like to distinguish several papers where the connection between Minkowski functionals and some functional characteristics of natural environs or liquid motion in them is demonstrated [21,22,23,24,25]. For example, in [25], Minkowski functionals were used to show that structure and water distribution impacts fungal colonization in soil. Another study investigated the intrinsic geometry change of oil reservoir voids during acid handling using topological methods and persistent homology theory [26].

The process of soil saturation with liquid leads to significant changes in soil properties. This process was discussed in the works [27,28,29]. In particular, [27] investigated the influence of the cycles of wetting and drying on the dynamics of the pore structure of three different soil types. In this study, the methodology for sample void intrinsic geometry analysis [30,31,32,33] was used to quantitatively describe soil changes in samples in dry conditions upon wetting. With this issue, the tomographic images of several samples from different horizons in dry conditions were obtained. Then, all samples were moistened and tomographed for the second time.

In the present investigation, a morphological “opening” operation was used to obtain a cumulative and non-cumulative pore size distribution of Minkowski functionals [24,33]. This method was chosen because it is necessary to obtain three-dimensional images of pores of various sizes to calculate the distributions of topological characteristics. Additionally, in this paper it is proposed to analyze voids Betti numbers in addition to Minkowski functionals due to the visible interpretation: zero Betti number is associated with the number of individual pores and the first Betti number is associated with number of tunnels in soil sample.

2. Materials and Methods

The Phaeozem albic (WRB 2014, version of 2015) from Vladimirskoe Opolie was the object of investigation. Vladimirskoe Opolie is situated in Russian valley northwest of Vladimir, on the west high bank of river Kliaz’ma, to the southeast of the Moscow glaciation boundary. The average long-term precipitations total is 575 mm and the evaporation capacity is more than 400 mm per year. Some physical properties of these soils are represented in

Table A1. Physical and chemical properties of silty loam Phaeozem albic.

Depth, cm	Horizons	Soil Texture Diameter, mm			Density, g/cm3	FC, % vol.	Sat. Hydraulic Conductivity, cm/Days
		<0.002	0.002–0.05	>0.05
0–5	Arable	17.39	80.66	1.95	1.1	37.3	60
0–10		17.35	80.21	2.44	1.16	37.2	58
10–20		17.21	80	2.79	1.21	37	52
20–30	A2	17.03	81.62	1.75	1.33	38.6	26
30–40		16	82.43	1.57	1.36	38.2	32
40–50	AB	17.35	81.76	0.89	1.33	37.4	35
50–60		17.32	82.09	0.59	1.39	35.3	35
80–100	B2	18.4	81.1	0.5	1.4	37.6	28.5

. The size of investigated objects was determined by the laser diffraction method, in which the laser particle analyzer Analysette 22 Comfort (FRITCH, Germany) was used. The bulk density was determined with sample rings, the filed capacity was determined in the field [34].

In Phaeozem albic the subsoil compaction appears at 20–30 cm depth in terms of consolidation at 20–40 cm depth. The consolidation graduates into horizons compacted by evolutionary soil processes in deeper layers. The compaction does not go beyond the critical values of bulk density for loamy cultivated soil (>1.4 g/cm³) [35]. However, the differences are essential compared to arable layers which have bulk density on optimum level. This occurrence stands to be mentioned as a characteristic one for clay-loam intensively used ploughing soils.

2.1. Sample Preparation and Wetting Process

The mold for monoliths was constructed using a medicine syringe with 20 mL volume. The syringe hub was taken out and the syringe was cut down to a height of 40 mm. For sample collection, soil was cut until it formed a rail of 10 cm depth around the monolithic sample. Thus, a “pocket” of soil was made in the form of truncated cone, where the upper base diameter was congruent to the syringe diameter and the lower base diameter was 2–3 cm wider than the syringe diameter. The form of the truncated cone provided stability and reduced the risk of monolithic cracking or breaking. After preparation procedures, the syringe was pushed on the monolith by undercutting the truncated cone to the cylindrical form every 5–10 cm of syringe moving. After the monolith was fully immersed into the syringe, the sample was cut in and evened from the bottom on the periphery of the syringe and was sealed with a sticky tape.

Thirteen monolithic samples were collected for laboratory investigation with X-ray computed tomography. Samples were collected from different horizons of Phaeozem albic soil. A visualization of voids from all samples in dry conditions can be found in Figure 1.

X-ray computed tomography was carried out under two different saturation conditions: air-dried and saturated beyond the field capacity. After the first air-dry tomography scanning, the sample was moistened by supplying excess water through the base standing on several filter papers. The imbibition continued for 7 full days with permanent watering until the water table rose up to the top base of the specimen. Then, the sample was free drained on a sandy base reaching the constant weight without evaporation. Then, a 3D image of the specimen was acquired through microtomography, with a time acquisition of 2–3 h. Soil humidity was constant during imaging.

2.2. X-ray Computed Tomography and Image Analysis

The tomographic investigation of pore space structure was carried out with a SkyScan 1172 (Belgium) microtomographer with different energy settings for dry (100 kV, 100 μA, Al + Cu filter) and wet (70 kV, 129 μA, Al filter) samples. The X-ray source was a Hamamatsu 100/250 with 10 W output. The imaging parameters were chosen so as to minimize costs (the running time of X-ray tunnel and CMOS device) of the microtomographer, but to obtain sufficient final image quality. Soil samples consisted of several vertical sections which were acquired in one hour.

The projections were reconstructed onto a 3D image with the Feldkamp algorithm [36] implemented in NReckon software. Voids visualization of each sample in dry and moisture-saturated conditions are provided in Supplementary Materials (Figures S1–S11). We obtained 3D sample images with 16-bit arrays and a resolution of 8 micrometers in each direction. After median filtering with a 3 by 3 kernel, pores were identified via Fiji [37] implementation of minimum cross entropy thresholding method [38], which showed the most consistent results. We will only consider the visible porosity because the pores smaller than several voxels in size are not resolved.

Two-dimensional tomographic images of one soil sample are presented in Figure 2. The sample was acquired in dry and moisture-saturated conditions. The changes in void structure are visible—some flat cracks are closed, and cylindrical shaped channels have narrowed.

3. Minkowski Functionals and Betti Numbers

3.1. Theory

Let $X$ be an object which is bounded by the smooth surface $\delta X$ in Euclidean d-dimensional space. To describe the geometric and topological properties of this object, the integral geometry enables us to determine $d+1$ topological invariants that are named Minkowski functionals [39,40]. Here, consideration is given to three-dimensional tomographic images of soil samples. Thus, only 3-dimensional space will be reviewed. The object of investigation was the pore space because its structure mainly determines the soil properties. In this case, the Minkowski functionals correspond to pores volume, surface area of pores, and integral mean curvature of pore–soil interface and Euler-Poincare characteristic of pore space:

$$M_0\left(X\right)=V\left(X\right),$$

(1)

$$M_1\left(X\right)={{\displaystyle \int }}_{\mathsf{\delta }X}dS=S\left(X\right),$$

(2)

$$M_2\left(X\right)={{\displaystyle \int }}_{\mathsf{\delta }X}\left(\frac{1}{r_1}+\frac{1}{r_2}\right)dS=C\left(X\right),$$

(3)

$$M_3\left(X\right)={{\displaystyle \int }}_{\mathsf{\delta }X}\frac{1}{r_1{r}_2}dS=2\mathsf{\pi }\mathsf{\chi }\left(\mathsf{\delta }X\right)=4\mathsf{\pi }\mathsf{\chi }\left(X\right),$$

(4)

where $r_1$ and $r_2$ are the main radii of surface curvature $dS$, while $\chi \left(\delta X\right)$ and $\chi \left(X\right)$ are the Euler-Poincare characteristics for the surface and convex body, respectively. To compare values of Minkowski functionals obtained from different samples, we used normalization. For example, if the size of the 3-dimensional image was 1000 × 1000 × 1000 voxels, then the values of functionals were divided by 10⁹.

The Euler-Poincare characteristic for the convex body is an integral estimate of topological complexity, and can be determined as an alternative summary of Betti numbers:

$$\chi \left(X\right)=b_0\left(X\right)-b_1\left(X\right)+b_2\left(X\right),$$

(5)

where Betti numbers can be interpreted in the context of soil pore space in the following way [31,41]: $b_0$ is the number of connected components (individual pores), $b_1$ is the number of “tunnels” (or redundant connections, loops, torus—all types of pores that cannot be compacted into a dot), and $b_2$ is the number of solid clusters that are completely surrounded by pores (supposed to be equivalent to 1).

Only two phases (pores and solid skeleton) are rationalized in the investigated sample, so, according to [42,43], Minkowski functionals and Betti numbers were calculated for the pore space. The water inside the specimens in moisture-saturated condition was not considered. Values for the solid phase can be determined using an analytic approach.

3.2. Algorithms

The 3-dimensional image of voids can be represented in terms of cubical complexes [44,45]. The additivity of Minkowski functionals affords the use of combinatorial algorithms to calculate them. There are several algorithms which enable the calculation of Minkowski functional values [40,46,47]. In [46], an algorithm based on calculating the number of different local configurations of 2 × 2 × 2 voxel was proposed, while in [40] an algorithm based on calculating the number of vertices, edges, faces and voxels in cubical complex was presented:

$$V\left(X\right)=n_c,$$

(6)

$$S\left(X\right)=-6n_c+2n_f,$$

(7)

$$2C\left(X\right)=3n_c-2n_f+n_e,$$

(8)

$$\mathsf{\chi }\left(X\right)=-n_c+n_f-n_e+n_{\mathsf{\nu }},$$

(9)

where $n_c$ is the number of cubes, $n_f$ is the number of faces, $n_e$ is the number of edges and $n_{\nu }$ is the number of vertices in the cubical complex which is associated with sample voids.

In this study, we used [47] software to compute Minkowski functionals. Zero and second Betti numbers (numbers of connected clusters of pore space and solid matrix) were calculated using basic functions of MATLAB [48] and ImageJ [49]. The first Betty number was calculated using the following equation:

$$b_1\left(X\right)=b_0\left(X\right)+b_2\left(X\right)-\mathsf{\chi }\left(X\right).$$

(10)

3.3. Mathematical Morphology: Erosion, Dilatation, and Opening Operations

The morphological opening operation used in this study relates to the type of methods which is called the mathematical morphology. Classical mathematical morphology was proposed in [50], where the main morphological operations of erosion and dilatation were introduced. Mathematical morphology was used for tomographic soil images [20,33,51] analysis. A detailed mathematical morphology methods review describing their applicability for investigations of 2D and 3D images of natural structures is provided in [52].

The erosion operation leads to object boundaries expansion. Conversely, the dilatation leads to boundaries narrowing (Figure A1 and Figure A2). An opening operation was obtained by successive usage of erosion and dilatation operations with structural unit identical form and size for objects:

$$X\circ B_d=\left(X\ominus B_d\right)\oplus B_d,$$

(11)

where $B_d$ is the structural spherical shaped unit with diameter $d$, and symbols $\ominus and \oplus$ define erosion and dilatation operations, respectively. As a result of the opening operation usage, objects which are smaller than the size of structural unit disappear. Thus, the dependences $M_i\left(d\right)$ and $b_i\left(d\right)$ represents the cumulative distribution of $M_i$ and $b_i$ over pore sizes. Non-cumulative distributions of Minkowski functionals and Betti numbers can be carried out as a derivative of the cumulative distributions.

In Figure 3 demonstrates the results of the morphology opening operation with increasing structural discus-shaped unit appliance to the sample pore space image. The generation of similar functionals distributions for the pore space of samples in dry and wet conditions enabled the quantitative description of change in pore space structure due to soil wetting. The advantage of this approach is that there is no need to separate the pore space into individual pore objects [5].

4. Results and Discussion

4.1. Topological Analysis of All Samples

The topological invariants (Euler-Poincare characteristic and Betti number) of all samples in our investigation are reported in

Table 1. Soil samples used in this study, indicating the occurrence depths and horizons. The three columns on the right represent the relationship between the Euler-Poincare characteristic and Betti numbers of the pore space under dry and wet conditions.

Sample	Depth, cm	Horizon	Topology	${\mathit{\chi }}^{\mathit{w}\mathit{e}\mathit{t}}\mathit{v}\mathit{s}.{\mathit{\chi }}^{\mathit{d}\mathit{r}\mathit{y}}$	${\mathit{b}}_0^{\mathit{w}\mathit{e}\mathit{t}}\mathit{v}\mathit{s}.{\mathit{b}}_0^{\mathit{d}\mathit{r}\mathit{y}}$	${\mathit{b}}_1^{\mathit{w}\mathit{e}\mathit{t}}\mathit{v}\mathit{s}.{\mathit{b}}_1^{\mathit{d}\mathit{r}\mathit{y}}$
1	0–10	Arable	Irregular	>	>	<
2	0–10	Arable	Normal	>	<	<
3	0–10	Arable	Irregular	>	>	<
4	10–20	Arable	Irregular	>	>	<
5	10–20	Arable	Irregular	>	>	<
6	20–30	A2	Normal	>	<	<
7	20–30	A2	Normal	<	<	<
8	30–40	A2	Normal	<	<	<
9	30–40	A2	Normal	<	<	<
10	40–50	AB	Normal	<	<	<
11	40–50	AB	Normal	-	-	-
12	80–90	B2	Normal	<	<	<
13	80–90	B2	Normal	<	<	<

. For samples from the arable horizon, an increase in the Euler-Poincare characteristic is observed during the wetting process. On the other hand, for samples from the A2, AB, and B2 horizons, the Euler-Poincare characteristic decreased during wetting (except for sample 6). All samples were divided into two groups based on mechanism of topological change of pore space during wetting:

• Samples with a “normal” topology change (samples from 20–90 cm depth, horizons A2, AB, and B2). Both Betti numbers (number of connected clusters and number of tunnels) decreased during swelling. The Euler-Poincare characteristic change was determined by the correlation between closed pores and closed tunnels: In the samples from the above-mentioned horizons the number of closed pores was bigger than the number of closed tunnels: ($b_0^{\mathit{dry}}-b_0^{\mathit{wet}}>b_1^{\mathit{dry}}-b_1^{\mathit{wet}}$), so the Euler-Poincare value decreased. These samples were taken from the horizons which had not been agriculturally exploited, so a denser (>1.3 g/cm³) and more stable structure was preserved. Upon wetting, some tunnels are closed, or at least partly closed, allowing them to be subsumed to the category of individual pores. In one of the samples taken from 20–30 cm depth from the A2 horizon (sample 6), the number of closed pores was smaller than the number of closed tunnels and the Euler-Poincare value increased. The correlation can be explained by the individual structural features of the samples, particularly by the large number of small tunnels which closed during wetting (Figure S17).
• Samples with an “irregular” changing topology (0–20 cm, arable horizon). In these samples, the $b_1$ (number of tunnels) decreased, but the $b_0$ (number of connected clusters) increased during the wetting process. In this case, an increase in the number of small pores was observed (this was also proved by the increasing integral mean curvature in that range). Apparently, during the swelling, some pores become smaller but not enough to close. Another explanation can be a mechanism whereby small pores do not close up because they are filled with an X-ray transparent substance, for example, clamped pendular water [53]. In samples from the arable horizon at a depth of 10–20 cm, a number of tunnels narrowed slightly, which can be explained by the presence of plant roots in most tunnels. For the majority of agricultural crops in a mild climate, about 50% of all plant roots are accumulated at a depth of 8–20 cm [54]. These roots can keep tunnels open during soil wetting. It should be noted that in one of the samples from the arable horizon at a depth of 0–10 cm (sample 2), the change in the topology of the pore space occurred in agreement with the “normal” change (Figure S13).

4.2. The Detailed Analysis of the Sample from the B2 Horizon

Here, we discuss the pore size distribution of Minkowski functionals and Betti numbers in the sample from the B2 horizon. The sample was taken from a depth of 80–90 cm and its pore space image is presented in Figure 4. The dependence graphs for all samples are shown in Supplementary Materials (Figures S12–S22). The dependencies analysis from Figure 5 suggests that, during soil wetting, the total sample porosity narrowed across the entire range of pore sizes; furthermore, the specific surface area of pores narrowed too. The integral mean curvature of pore spaces narrowed by less than 0.4 mm, which can be explained by the small number of pores closing. These small pores have a large specific surface curvature. The large pores and tunnels reduced slightly in size and did not change their form, so their curvature changed minimally.

More detailed changes which take place in the void’s structure are illustrated by the distribution of the Euler-Poincare characteristic and Betti numbers. The larger the Euler-Poincare value is, the stronger individual pores dominate over tunnels in a void’s structure. The structures where tunnels prevail over individual pores have negative values of the Euler-Poincare characteristic.

In a case sample from B2 horizon, a Euler-Poincare characteristic of voids narrows during soil swelling across the entire range of pore sizes. This means the voids become more bound together topologically. This can be proved by Betti numbers analysis, which enables visual interpretation. A zero Betti number is associated with the void’s connected clusters, whereas the first Betti number is associated with numbers of tunnels in the solid phase. The individual pore spaces ($b_0$) are interstices that do not have exits past the sample boundaries. These pores are closed. The tunnels ($b_1$) are open pores with more than one exit past the sample boundaries (perforated pore). It can also be a closed pore with a form topologically identical to torus.

The distribution analysis suggests that the dry sample contains more individual pores or tunnels than the wet sample, with a difference in the range of 0–0.4 mm. Nevertheless, pores are more likely to become closed than tunnels: $b_0^{\mathit{dry}}-b_0^{\mathit{wet}}>b_1^{\mathit{dry}}-b_1^{\mathit{wet}}$. It is worth noting that the aforementioned correlation is not universe for all samples. The second Betti number equals 1 (it is associated with the connected clusters of soils phase which are surrounded by voids) and is not represented in graphs.

5. Conclusions

Here we presented the methodological base of the void structure quantitative description. The methodology is based on the methods from integral geometry, topology, and mathematical morphology. Methodology can be applied to media with complex internal structures at different scales, including but not limited to soils, sedimentary rocks, foams, ceramics, and composite materials.

The void structure’s evolution of 13 monolithic samples during soil wetting was quantitatively illustrated. To quantitatively demonstrate the changes in the structure of the pore space, an analysis of the evolution of the Minkowski functionals and Betti numbers was carried out. The results suggest that quantitative topological analysis can be used to describe either the void structure or the changes in internal structure of voids which are initiated by some external effect or process. We found that the topology of the void space in the process of soil saturation can change in two ways. A “normal” change in topology implies a decrease in the zero and first Betty numbers of the pore space, and is observed in soils from the A2, AB, and B2 horizons at a depth of 20–90 cm. An “irregular” change in the pore space was found for soils from the arable horizon at a depth of 0–20 cm. The zero Betty number of the pore space for these samples increased and then decreased during wetting.

Eventually, we plan to transition from quantitative characteristics of pore space to the functional characteristics of hydro-physical properties of soil. The future target of research is the degraded soils. The possibility of association of Minkowski functionals and Betti numbers with physical properties such as field capacity, wilting moisture content, bubbling pressure, hydraulic conductivity, and other physical pedological characteristics will be investigated.