Multifractal Analysis of 3D Correlated Nanoporous Networks

Abstract

In this study, we utilize Monte Carlo methods and the Dual Site-Bond Model (DSBM) to simulate 3D nanoporous networks with various degrees of correlation. The construction procedure is robust, involving a random exchange of sites and bonds until the most probable configuration (equilibrium) is reached. The resulting networks demonstrate different levels of heterogeneity in the spatial organization of sites and bonds. We then embark on a comprehensive multifractal analysis of these networks, providing a thorough characterization of the effect of the exchanges of nanoporous elements and the correlation of pore sizes on the topology of the porous networks. Our findings present compelling evidence of changes in the multifractality of these nanoporous networks when they display different levels of correlation in the site and bond sizes.

1. Introduction

Natural or synthetic porous materials are characterized by their unique feature of having hollow spaces or pores immersed in a solid or semi-solid matrix [1]. These pores act as reservoirs for wetting fluids (air, water, nitrogen) and non-wetting fluids (petroleum, mercury) [2]. The diverse physicochemical phenomena occurring within these hollow spaces are determined by the characteristics of the pores, such as size, shape, or geometry; interconnection; quantity and location; or spatial distribution. These unique characteristics, in conjunction with the properties of the solid phase and the fluid dynamics within the hollow space, confer upon porous materials their distinct properties, such as mass, density, porosity, the volume of the pore space, pore surface area, roughness, permeability, and tortuosity [3,4].

The porosity of a material is not always evident to the naked eye. In these cases, techniques such as electron microscopy, nitrogen adsorption, helium pycnometer, mercury porosimetry, X-ray diffraction, and nuclear magnetic resonance techniques [5,6,7,8] help to explore the hollow space of materials. Organizations such as the International Union of Pure and Applied Chemistry (IUPAC) [9,10] have proposed standardizing symbols and terms and standardizing both characterization methodologies and how to present the information obtained. In particular, IUPAC proposed the nitrogen adsorption–desorption protocol to characterize porous materials. This protocol is essential for equipment manufacturers, as it helps them implement the software (based on the theory associated with the description of some textural properties of porous solid materials), which allows them to calculate values such as pore size, pore size distribution function, pore volume, and surface area.

Many real porous media are heterogeneous, with chemical impurities and crystalline imperfections being significant causes of this heterogeneity. These factors lead to variations in properties such as connectivity, symmetry, pore size, energy potentials, and vapor–solid adsorption [1]. The irregular manifestation of these properties from one point to another in the porous medium underscores the system’s complexity [3]. From the above, it is deduced that the complexity of porous media has promoted the development of tools that allow for investigation of each of their characteristics from different approaches. A crude and simple way to represent porous media is through models. To do this, the initial and oldest models considered the hollow space of a material as a set or bundle of independent cylinders (or spheres) of various sizes that were not interconnected. Later, models of one-dimensional, two-dimensional and three-dimensional networks were proposed. These networks, often referred to as “lattices”, consist of interconnected “sites” (also called “nodes”, “junctions”, “cavities”, and “hollow spheres”) and “bonds” (also called “capillaries”, “channels”, “throats”, “necks”, “windows”, “hollow cylinders”, and “links”) were proposed. For instance, 1D networks are characterized by the fact that the diameter along their axis changes in an irregular manner (they are a linear chain of alternately connected sites and bonds). In 2D networks, the difference lies in the way in which their porous entities are connected (connectivity (C) is the mean number of bonds surrounding a site [11] and measures the degree to which a porous structure is interconnected [4]); they can be of the types tri-tetragonal ($C=2.66$), hexagonal ($C=3$), hexa-triangular ($C=6$), di-tetragonal ($C=3$), triangular ($C=6$), square ($C=4$), tetra-triangular ($C=6.6$), Bethe ($C=3$) [12,13], honeycomb ($C=3$), and Kagomé network ($C=4$) [14]. In the examples above, the networks all share a common feature—symmetry. For instance, the hexagonal, di-tetragonal, honeycomb, and Bethe networks all have a connectivity of three. However, the differences in their topology, as well as the spatial arrangement of sites and bonds, give them unique properties. It is important to note that, while networks can be irregular, regular networks are more straightforward to study [1]. Finally, 3D networks are considered a set of hollow sites placed at the nodes of a lattice interconnected through hollow cylindrical bonds.

As can be seen, porous networks are characterized by having intrinsic parameters, such as dimensionality, connectivity, and topology. In general, the proposed models consider the porous medium as a network of large spaces connected to smaller conduits of specific geometries, thus reducing the degrees of freedom by discriminating many fundamental physicochemical properties [4]. However, one of the benefits of modeling is that “tailored models” can be proposed according to the application we wish to explore. For example, Rodríguez et al. [15] used log-normal size distribution functions of sites (cylindrical disks) and bonds (rectangular prisms) to construct a 2D correlated network. In it, they simulated the growth of bubbles. In turn, Rojas et al. [16] used Gaussian twofold distributions (spherical sites and cylindrical bonds) to construct three-dimensional cubic porous networks (variable connectivity, C = 2 to 6). The phenomenon of nitrogen adsorption and desorption was simulated on these heterogeneous networks.

On the other hand, fractal and multifractal analysis are valuable tools for describing irregular shapes without a characteristic scale. Multifractal analysis is appropriate and highly practical for image analysis. It can effectively capture and characterize the complex scaling behaviors and local variations in inhomogeneous images. More specifically, estimating local exponents or generalized fractal dimensions makes it possible to estimate the complexity of images. In other words, with this information, it is possible to characterize the inhomogeneous scaling properties of images, which are essential for understanding the local and global irregularities of the image. In this way, numerous studies have used these analysis techniques to describe various porous media and the processes that occur inside them [17,18,19,20]. For example, the microstructure of bubbles inside gaseous porous sediments has been characterized [17]; scanning electron microscopy images of human dentin have been analyzed to find its porous structure by multifractal analysis [21]; and different fluid displacement patterns inside 3D porous media have been found with information from techniques such as mercury porosimetry, nuclear magnetic resonance, electron microscopy, and X-ray computed tomography [18].

In this work, we use the Monte Carlo [22] method and the Dual Site-Bond Model to simulate constructing three 3D nanoporous networks with different degrees of correlation. Random exchanges of their elements (sites and bonds) were performed in each 3D network, which was built under the assumptions of the DSBM model. Throughout the simulation, the structure of the networks changed (evolved) through different configurations until the most probable (equilibrium) configuration was reached. In this way, each specific number of network element exchanges (Monte Carlo steps) was recorded, obtaining a series of 12 intermediate networks for each configuration initially established. Then, a multifractal analysis of the simulated network series was performed to evaluate the spatial organization of these porous networks. We find that networks with moderate correlation tend to exhibit greater changes in multifractality, while the higher the correlation, the more stable the multifractality as the networks evolve.

The structure of the paper is as follows. In Section 2, we present the main aspects of the methodology employed to simulate the networks, as well as a brief explanation of the multifractal method used to characterize the networks. Next, in Section 3, we present the results of different levels of multifractality observed in the simulated data. Finally, we include our conclusions in Section 4.

2. Methodology

2.1. Models

When describing a complex and intricate system, the most straightforward approach is to propose a simple ideal model. This model should allow us to determine the most significant capabilities of the system using a minimal number of basic parameters. In the case of porous media, the model needs to be broad enough to encompass all the specifics of a wide variety of porous materials. The basic parameters, such as the size distributions of the hollow elements and connectivity, are crucial. They serve as the foundation from which all other properties can be derived through the moments of these distributions (first moment: mean radius; second moment: surface; third moment: porosity; fourth moment: permeability) [23]. The nature of the porous medium largely determines the choice of pore size distribution shape; it can be uniform, normal, inverse-power, or Gaussian [4,24,25]. The key concept to bear in mind is the idea of duality; that is, accessing the network’s large hollows (sites) must be done from the small hollow passages (bonds).

A simple alternative model that describes the nature of porous media is the so-called “Dual Site-Bond Model (DSBM)”, which was proposed by Mayagoitia et al. [26,27]. The DSBM is a discrete model that considers the spatial correlation among pore sizes, thus allowing for the construction of pore networks with different structures. It is possible to visualize the network as a collection of spheres (sites) and cylinders (bonds) of sizes $R_S$ and $R_B$, respectively, that are also alternately connected [28], and whose collections are represented by their respective size distribution functions, $F_S\left(R\right)$ and $F_B\left(R\right)$, which must be normalized and can be expressed in terms of the total number of porous elements of each class: sites or bonds (Figure 1). Therefore, the probabilities of finding a site or bond of size less than or equal to R can be represented as $S\left(R\right)$ and $B\left(R\right)$, respectively:

$$S\left(R\right)={\int }_0^R{F}_S\left(R\right)dR\mathrm{and}B\left(R\right)={\int }_0^R{F}_B\left(R\right)dR.$$

(1)

The DSBM also establishes that if we want to build a self-consistent porous network, it is necessary to comply with the following Construction Principle (CP): “The size of a site is greater than or equal to the size of the bonds to which it is connected or vice versa, a given bond must have a size less than or equal to the size of the sites between which it is located” [29,30]. Based on this principle, it is possible to state two laws that guarantee that the CP is satisfied throughout the porous network:

$$B\left(R\right)-S\left(R\right)\ge 0,\mathrm{for}\mathrm{all}R.$$

(2)

This first law (Equation (2)) implies that the distribution function of the sites must always be located on the right or at most completely overlap with that of the bonds ($b1\le s1$ and $b2\le s2$; see Figure 1). Even if there is an appropriate collection of sites and bonds, that is, the first law is fulfilled both in size and in number, the CP can be violated if there is a certain degree of overlap $\Omega$ (common area) between the two distributions (Figure 1d,f), since there would be some bonds larger than some sites that, when incorrectly assigned (connected to smaller sites), would automatically violate the CP. To avoid this inconsistency, Mayagoitia et al. [26,27] proposed a correlation function $\varphi (R_S,R_B)$ that represents a topological correlation function between the sizes of sites and interconnected bonds and implies that the event of finding a size $R_S$ for a site connected to an $R_B$ size bond is not independent. In this way, the joint probability density, $F(R_S,R_B)$, of finding a site with size $R_S\in (R_S,R_S+d{R}_S)$ connected to a bond with size $R_B\in (R_B,R_B+d{R}_B)$ [31,32] is:

$$F(R_S,R_B)=F_S\left(R_S\right)F_B\left(R_B\right)\varphi (R_S,R_B),$$

(3)

where the correlation function $\varphi$ is introduced; if $\varphi (R_S,R_B)=1$, the previous events would be independent for all $R_S$ and $R_B$. Then, the network would be built entirely randomly; however, if $\varphi (R_S,R_B)\ne 1$, the events would be correlated. The presence of correlation leads to the second law:

$$F(R_S,R_B)=0\mathrm{for}{R}_S<R_B.$$

(4)

The second law (Equation (4)) implies that $\varphi (R_S,R_B)=0$ is local and expresses the fact that “the size of a site is greater than or equal to the size of the bonds to which it is connected” [28]. The most plausible form of $\varphi (R_S,R_B)$, that is, the one that allows the most significant degree of randomness in the construction of the porous network without violating the CP, is the following [26]:

$$\varphi (R_B,R_S)={\displaystyle \frac{exp\left(-\underset{B\left(R_B\right)}{\overset{B\left(R_S\right)}{\int }}\frac{dB}{B-S}\right)}{B\left(R_B\right)-S\left(R_B\right)}}={\displaystyle \frac{exp\left(-\underset{S\left(R_B\right)}{\overset{S\left(R_S\right)}{\int }}\frac{dS}{B-S}\right)}{B\left(R_S\right)-S\left(R_S\right)}}=\varphi (R_S,R_B).$$

(5)

Equations (5) and (3) are now well-determined and illustrate that the degree of overlap between $F_S\left(R\right)$ and $F_B\left(R\right)$ is related to the correlation function (correlation is understood here as how sites and bonds are grouped or connected along the porous network as the overlap of the size distribution functions changes without violating the CP). Figure 1a schematizes a porous medium whose site and bond sizes are well-differentiated by the lack of overlap ($\Omega =0$) between the distribution functions (Figure 1b). Using those sizes in network construction could only result in an uncorrelated network. In contrast, Figure 1e shows the extreme case, where the bonds are almost the same size as the sites, so the overlap is almost complete ($\Omega \approx 1$) between the size distribution functions (Figure 1f). Using those sizes in construction could only result in a highly correlated network. Finally, Figure 1c shows the intermediate case where the size of the sites and bonds is variable (from very small to very similar), so the size distribution functions (Figure 1d) partially overlap ($0<\Omega <1$). Using those sizes in network construction could only generate a moderately correlated network. It is important to note that the degree of correlation in the porous media is a complex and decisive factor that influences many capillary processes. For example, a higher correlation can lead to a more rapid invasion of wetting or non-wetting fluids. In contrast, a lower correlation favors delayed invasion or entrapment of fluids in the porous medium [11,33].

2.2. 3D Correlated Nanoporous Network Simulations

Despite the prevalence of heterogeneity in real porous media, we have made a deliberate decision in this work to limit the study variables and model the porous medium in the simplest way possible. This approach, which constructs the porous medium as a regular network while maintaining the dual idea of porous elements and regularity in connectivity, allows us to focus our research and ensures rigor in our analysis. The 3D correlated nanoporous networks with constant connectivity ($C=6$ in this work) were constructed using the premises of the Dual Site-Bond Model (DSBM). For this, normalized Gaussian functions were used to represent the pore-size distribution of sites ($F_S\left(R\right)$) and bonds ($F_B\left(R\right)$). The nanopore-size ranges of $F_S\left(R\right)$ and $F_B\left(R\right)$ span from ${\overline{R}}_S-3\sigma$ to ${\overline{R}}_S+3\sigma$ and ${\overline{R}}_B-3\sigma$ to ${\overline{R}}_B+3\sigma$, respectively, where ${\overline{R}}_S$ (${\overline{R}}_B$) represents the mean value of the sites (bonds), and $\sigma$ is the standard deviation ($\sigma =1.5$ nm in this work).

Table 1. Parameters of the two-fold Gaussian distributions used in constructing 3D correlated nanoporous networks using the DSBM. $\Omega$ is the overlap between $F_B\left(R\right)$ and $F_S\left(R\right)$.

Network	${\overline{\mathit{R}}}_{\mathit{B}}$ (nm)	${\overline{\mathit{R}}}_{\mathit{S}}$ (nm)	$\mathbf{\Omega }$
$NA$	22.5	24.0	0.62
$NB$	22.5	23.5	0.74
$NC$	22.5	23.2	0.82

shows the parameters of the Gaussian distributions used for the in silico construction of the $NA$, $NB$, and $NC$ nanoporous networks.

The simulated networks can be visualized as cubic lattices of $L^3$ sites and $C/2\times L^3$ bonds ($L=100$ in this work). The sites were placed on the lattice nodes, and the bonds were placed between them such that a site had six bonds surrounding it. Periodic boundary conditions were considered for each of the external faces of the cubic network.

The algorithm used to simulate 3D nanoporous networks is the following [11]:

• From the collection of sites of size R represented by $F_S\left(R\right)$, $L^3$ sites were randomly sampled and placed on the nodes of the cubic lattice.
• From the collection of bonds of size R represented by $F_B\left(R\right)$, $C/2\times L^3$ bonds were randomly sampled and placed on the nodes of the cubic lattice.
• Inconsistencies or violations to the CP of the DSBM were accounted for in the initial cubic lattice made up of the sites and bonds and were corrected through a random exchange of its elements; that is, two sites were chosen at random from the network, and an attempt was made to exchange from their nodal positions. The exchange was carried out only when both sites were greater than or equal to all the bonds that surrounded them (the CP of the DSBM is fulfilled). The same was done with the bonds, except that the bonds had to be smaller than the sites surrounding them.
• Once there were zero CP violations, a series of network transitions (“shake the network”) were performed as many times as necessary, going through different possible network configurations, until the most probable configuration (equilibrium) was reached. The network was “relaxed” by randomly exchanging its elements (the exchange of porous entities is a crucial process, carried out only when the CP is not violated in each exchange attempt). For this, a Monte Carlo Step (MCS) was defined as $L^3+C/2\times L^3$ exchange attempts of network elements (sites and bonds). Each specific number of network element exchanges (MCSs) was recorded, obtaining a series of 12 subsequent networks ($NA-i$, $NB-i$, and $NC-i$) for each initially established configuration ($NA$, $NB$, and $NC$). The letter i is an index representing a specific network configuration generated according to MCSs, i.e., $i\in \{1,2,3,...,12,13\}$.

  Table 2. Number of attempts to exchange nanoporous elements (MCSs) to achieve different configurations of the $NA$, $NB$, and $NC$ networks. The random exchange of elements transforms these networks into the specific configurations $NA-i$, $NB-i$, and $NC-i$, respectively. i is the letter used to identify the network number within a set of networks with the same ${\overline{R}}_B$ and ${\overline{R}}_S$ values.

  i	MCS	i	MCS	i	MCS
  1	0	6	60,000	11	300,000
  2	5000	7	80,000	12	400,000
  3	10,000	8	100,000	13	500,000
  4	20,000	9	150,000
  5	40,000	10	200,000

  presents the MCS values used to generate different network configurations as randomly as possible. In this study, the determination of sizes of the nanoporous elements (sites and bonds), the selection of the nanoporous elements, and the selection of the nodal positions of the lattice for exchanges between the elements were all carried out randomly. This was achieved by generating an algorithm of pseudo-random numbers [11,22].

2.3. Two-Dimensional MF-DFA (2D MF-DFA)

The following five steps describe the two-dimensional version of the multifractal detrended fluctuation analysis (MF-DFA) [34,35].

Step 1. Consider a self-similar surface represented by X(i,j), with $i=1,2,3,...,M$ and $j=1,2,3,...,N$. Let be $X_{m,n}$ the partition of the surface into $M_s\times N_s$ disjoint square segments of size $s\times s$, where $M_s=\lfloor M/s\rfloor$ and $N_s=\lfloor N/s\rfloor$. Each segment is denoted by $X_{m,n}=X_{m,n}(i,j)$ with $X_{m,n}(i,j)=X(u+i,v+j)$ for $1\le i,j\le s$, where $u=(m-1)s$ and $v=(n-1)s$. Given that M and N are not multiples of the length s, we ignore X’s last parts (at the end of the right and bottom).

Step 2. For each segment $X_{m,n}$, the cumulative sum $u_{m,n}(i,j)$ is defined as:

$$u_{m,n}(i,j)=\sum _{k_1=1}^i\sum _{k_2=1}^j\left\{X_{m,n}(k_1,k_2)-\overline{X}\right\},$$

(6)

where $\overline{X}$ is the average of the entire surface $X(i,j)$. This step is included to avoid a monotonic behavior in the cumulative sum.

Step 3. For each sub-surface $u_{m,n}$, we compute a local trend ${\tilde{u}}_{m,n}$ by fitting it with a linear bi-variate polynomial function:

$${\tilde{u}}_{m,n}(i,j)=ai+bj+c,$$

(7)

where a, b, and c are parameters that can be determined. Generally, it is possible to consider higher-order polynomials for the fitting function [36]. Next, the residual matrix is defined as follows:

$${ϵ}_{m,n}(i,j)=u_{m,n}(i,j)-{\tilde{u}}_{m,n}(i,j).$$

(8)

The detrended fluctuation function $F(m,n,s)$ of the segment $X_{m,n}$ is:

$$F^2(m,n,s)={\displaystyle \frac{1}{s^2}}\sum _{i=1}^s\sum _{j=1}^s{ϵ}_{m,n}^2(i,j).$$

(9)

Step 4. When q is non-zero, the generalized detrended fluctuation function of order q, denoted by $F_q\left(s\right)$, is computed as follows:

$$F_q\left(s\right)={\left\{{\displaystyle \frac{1}{M_s{N}_s}}\sum _{m=1}^{M_s}\sum _{n=1}^{N_s}{\left[F(m,n,s)\right]}^q\right\}}^{1/q}.$$

(10)

For the case $q=0$, the generalized detrended fluctuation function is [34]:

$$F_0\left(s\right)=exp\left\{{\displaystyle \frac{1}{M_s{N}_s}}\sum _{m=1}^{M_s}\sum _{n=1}^{N_s}ln\left[F(m,n,s)\right]\right\}.$$

(11)

Step 5. Finally, for multifractal surfaces, $F_q\left(s\right)$ follows a scaling power law:

$$F_q\left(s\right)\sim s^{h\left(q\right)},$$

(12)

where $h\left(q\right)$ is the Hölder exponent [37]. In multifractal formalism [38], the partition function is responsible for capturing the (multi)fractal properties by using the mass function $\tau \left(q\right)$, which is a nonlinear function of q [38,39]. For each q-value, we can obtain the $\tau \left(q\right)$ function through the following expression [34]:

$$\tau \left(q\right)=qh\left(q\right)-D_f,$$

(13)

where $D_f$ is the fractal dimension of the geometric support of the multifractal measure ($D_f=2$ in this work) [34,38]. From $\tau \left(q\right)$ (Equation (13)), we can obtain the generalized dimension $D_q$, the singularity strength function $\alpha \left(q\right)$, and the multifractal spectrum $f\left(\alpha \right)$ via the Legendre transform [34,37,38,39]:

$$D_q=\tau \left(q\right)/(q-1),$$

(14)

$$\alpha =d\tau /dq,\mathrm{and}f\left(\alpha \right)=q\alpha -\tau \left(q\right).$$

(15)

We consider the following quantities to characterize the multifractality:

• Width of the spectrum: distance between the maximum and minimum singularity strength:

  $$\Delta \alpha ={\alpha }_{max}-{\alpha }_{min},$$

  (16)

  where $f\left({\alpha }_{max}\right)\approx f\left({\alpha }_{min}\right)$.
• Main singularity strength ${\alpha }^{\ast }$:

  $$f\left({\alpha }^{\ast }\right)=max\left\{f\left(\alpha \right)\right\}.$$

  (17)
• The support, information, and correlation dimensions:

  $$D_0,D_1,\mathrm{and}{D}_2.$$

  (18)

3. Results and Discussions

We constructed three unique 3D nanoporous networks ($NA$, $NB$, and $NC$) with different levels of correlation using $L^3$ sites, $C/2\times L^3$ bonds (L = 100), and the parameters from Table 1. The resulting normalized size distribution functions, $F_S\left(R\right)$ and $F_B\left(R\right)$, are shown in Figure 2. The only difference between the three networks was the overlap ($\Omega$) between the size distribution functions (Figure 2). Figure 3 shows three snapshots of each network for different stages in the evolution. We observe that the networks exhibit topological changes as the number of MCSs increases. Each topology state is essentially due to the number of exchanges carried out on the network elements, i.e., the number of MCSs (see Table 2). In particular, the $NA$ network, constructed with low overlap between site and bond distributions, exhibits significant topological changes as they evolve (Figure 3a1–a3). Instead, we observe that the $NB$ and $NC$ networks undergo slight changes in their topological structure as they evolve (see panels b1–b3 and c1–c3 of Figure 3).

This higher correlation favors a “size segregation effect”, whereby larger sites and bonds form regions of larger elements. In contrast, smaller elements come together to form regions of small entities. As explained above, this effect becomes more critical as $\Omega$ increases [40]. It is evident that the role of overlap (${\Omega }_{NC}$ > ${\Omega }_{NB}$ > ${\Omega }_{NA}$) is pivotal, as it directly influences the size of the “size segregation patches”, making them more significant. It can also be noted in Figure 3 that the greater the number of exchanges of porous elements, the more the regrouping of porous entities of similar sizes is favored.

The multifractal spectrum, a tool that provides crucial information about pore size distribution and spatial arrangement [41], is a fundamental component of our research. We employ the 2D MF-DFA method to characterize the heterogeneity observed in the 3D nanoporous networks. This method, which involves analyzing slices of the network in two dimensions, allows us to capture the complex spatial distribution of the pores. Figure 4 illustrates our procedure for the network $NA-1$ case. Our analysis considers slices perpendicular to the three Cartesian axes, generating 100 images per axis and 300 for each network (see panels a and b1–b3 in Figure 4). We then constructed each image’s fluctuation function $F\left(s\right)$ across different scales s. The scaling behaviors, characterized by the local exponent $h\left(q\right)$ for different q values, are shown in panels c1–c3 of Figure 4. From this information, the principal representative magnitudes of the multifractal analysis were calculated, and their mean values (on the images or slices obtained from the three Cartesian axes) were considered for the characterization of the 3D networks.

In Figure 5, we show the results of the multifractal analysis for all the networks described in Table 1 and Table 2. We found that the multifractal spectra exhibit a wide range of values of local exponents $\alpha$, for which $f\left(\alpha \right)$ displays different levels (convex shape) with a clearly defined maximum. The function $f\left(\alpha \right)$ ranges from 0 to 2, which is to be expected since we are analyzing two-dimensional objects.

For NA networks (Figure 5a), we observe changes in the multifractal spectrum as the MCS increases, i.e., the range of scaling exponents (${\alpha }_{min}<\alpha <{\alpha }_{max}$) for which a fractal dimension $f\left(\alpha \right)$ exists shifts to the right, suggesting that the structure becomes more regular. Notably, for the network index $i>5$, the spectra tend to be increasingly overlapping (similar range), indicating that the distribution of the porous elements shows a saturation level from this value onward. The above is confirmed by the behavior of the local Hurst exponent ($h\left(q\right)$) in Figure 5b, where it is observed to be a monotonic decreasing function of q, and its value increases as the simulation evolves. In contrast, for the $NB$ networks, a slight shift of the spectrum to the right is identified as the network index increases (Figure 5c,d), while for $NC$ networks, no changes in the spectra are identified as the number of permutations increases (Figure 5e,f). Our results show that, for the set of simulated networks with the lowest overlap between the distributions ($\Omega =0.62$), the local exponents $\alpha$ are sensitive to the changes originated by the random exchanges of the elements. In contrast, the fractal dimensions associated with these local exponents display a similar behavior (almost identical convex shape). On the other hand, when the overlap increases ($\Omega =0.82$), the local exponent does not show sensitivity to increasing element exchanges, even when the number of MCSs is large.

To strengthen our analysis of the multifractality associated with nanoporous networks, we analyze the behavior of the generalized dimensions for three cases: the fractal or capacity dimension $D_0$, the information dimension $D_1$, and the correlation dimension $D_2$. Figure 6 shows the behavior of these dimensions regarding the network index for the three series of networks with different overlapping $\Omega$ values. For $D_0$, the three series of networks exhibit the same value ($D_0\approx 2$) for the space-filling capacity and are independent of the time evolution or the number of Monte Carlo Steps (MCSs). The evolution of $D_1$ is depicted in Figure 6b for the three series of networks. For the networks with the lowest overlap ($NA$), the dimension of the information grows from a value close to 1 to saturate at a value close to $1.6$ when the number of MCSs is large. On the other hand, both networks ($NB$ and $NC$) maintain a constant value ($D_1^{NB}\approx 1.7$, $D_1^{NC}\approx 1.8$) in the information dimension, indicating that the level of uniformity in the distribution of sites and bonds remained almost constant. It should be noted that, throughout the evolution, the highest value of the information dimension corresponds to the $NC$ networks, followed by the $NB$ and $NA$ networks. A similar behavior is observed for $D_2$, but with lower values in the three series of networks and a greater separation between their values along the evolution (see Figure 6c). These results for the correlation dimension mean that the dimensionality of the space occupied by sites and bonds undergoes essential changes when the networks come from configurations with lower overlap ($\Omega =0.62$) compared to cases where the overlap is higher ($\Omega =0.82$). Also, we remark that the generalized dimensions follow a clear descending order so that $D_0>D_1>D_2$ is satisfied, which confirms the multifractal character of the networks and their evolutions.

Figure 7a provides insight into the evolution of the singularity strength ${\alpha }^{\ast }$, i.e., the exponent with the most significant presence in the network. This exponent exhibits behavior similar to that observed for the generalized dimensions (see Figure 6b,c), revealing that networks constructed with the lowest overlap ($NA$) are the most sensitive to random element exchanges or MCSs. It should be noted that, throughout the simulation, ${\alpha }_{NA}^{\ast }<{\alpha }_{NB}^{\ast }<{\alpha }_{NC}^{\ast }$ is valid, which confirms that networks with high overlapping distributions lead to configurations with a higher level of regularity compared to the case where the correlation is lower. Regarding the width of the $\Delta \alpha$ spectrum (Figure 7b), we observe that the $NB$ networks exhibit the largest width throughout almost the entire evolution, followed by the $NA$ and the $NC$, with the width for the $NC$ being slightly smaller. However, at the end of the evolution ($i=13$), the $NA$ network reaches a greater width than those corresponding to $NB$ and $NC$.

4. Conclusions

In this manuscript, we have presented the in silico construction of 3D correlated nanoporous networks using Monte Carlo methods and under the premises of the DSBM. Our simulations of these networks have shown that the random exchange of porous elements (sites and bonds) causes the following effects on network topology: (i) in low-correlated networks, significant changes in the pore structure occur in the initial stages of the exchange; i.e., the dynamics of exchanges tends to form small clusters of porous elements of similar sizes (increases persistence); (ii) in highly correlated networks, from their initial construction, they are highly structured and are formed with massive clusters composed of elements of similar sizes. Thus, the exchange of their elements does not affect the topology, since the porous elements are always found in domains of similar sizes during their exchange. On the other hand, the multifractal analysis has shown that the conditions imposed by the level of overlap (or correlation) between the size distribution functions of the sites and the bonds are decisive for reaching configurations where different levels of heterogeneity are observed. When the correlation is relatively high, the multifractal spectra remain unchanged throughout the evolution, indicating the network is insensitive to element exchanges, whereas when the correlation is lower, the spectra exhibit shifts toward more persistence, showing heterogeneity distinct from the high-correlation case. Unlike other techniques for characterizing the porous structure of materials, multifractal analysis has made it possible to characterize nanoporous networks at different spatial scales. Multifractal analysis can identify the spatial variations inherent to the structures employing local densities and can concentrate the results of the behaviors in the multifractal spectrum. These results provide quantification of the presence of local and global irregularities (variations) present in a non-homogeneous porous structure with different levels of correlation, which is not easy to infer by traditional physical methods such as nitrogen adsorption [9,10], scanning electron microscopy [42,43], or mercury porosimetry [44].

The results presented here align with previous studies on the levels of multifractality observed in porous systems [41,45,46,47,48] and could improve our understanding of capillary processes [49] and fluid flow in porous media [20,48]. For example, in the Hg intrusion–retraction process, the shape of its curves could be expected to be similar in highly-correlated networks ($NC$). However, in low-correlated networks ($NA$), despite having the same size distribution functions, the spatial distribution of the pores is very different. This different topology could significantly impact the Hg invasion–retraction processes and the amount of Hg entrapped in the porous structure [11]. Our methodology, which includes 2D MF-DFA, a robust statistic that ensures more detailed characterization, has allowed us to examine pore networks from a multifractal perspective, and we have glimpsed crucial features and changes in spatial organization.