Suitability evaluation of an equivalent porous medium method for characterizing 2-D flow based on the fracture fractal dimension

Fractured rock mass, which is widely distributed in nature, is an important medium for fluid and solute transport. Fracture water is a significant source of water supply from a fractured geological medium (FGM) in many shallow buried bedrock distributed areas. Furthermore, with the increasing large-scale underground engineering activities, many countries pay more attention to the problems related to fracture water, such as groundwater resource extraction (Castellazzi et al. 2021), high-level radioactive nuclear waste disposal (Onoe et al. 2021), geothermal fluid exploitation (Guo et al. 2022), and underground crude oil storage facility evaluation (Li et al. 2014). The seepage model in the FGM can more effectively help to solve various practical problems in geological engineering analysis (Dong et al. 2019). Traditionally, the models are divided into two main categories: equivalent continuous medium model and discontinuous medium model, and the typical two models are an equivalent porous medium (EPM) and a discrete fracture network (DFN) for the two categories, respectively (Hadgu et al. 2017).

These two methods have different assumptions and applicable conditions and have their advantages and disadvantages. The EPM method is the most widely used for solving groundwater yield calculations at a field-to-regional scale. It considers that the permeability characteristics of fractured rock mass are similar to those of the porous continuous medium, and assumes that groundwater is distributed in all parts of rock mass (Wei et al. 2021; Pandey & Sharma 2023). The method does not need to know the specific distribution of the fracture network in the rock mass but only needs to know the statistical value of the equivalent permeability coefficient of rock mass (Scanlon et al. 2003). Therefore, the EPM model is a good choice for some seepage problems with severe weathering and fully developed joints in FGM areas (Zareidarmiyan et al. 2021). However, there are two main problems when using the EPM model. One is how to determine the suitability of the EPM model, and another is how to determine the equivalent permeability tensor of the FGM (Zeng et al. 2021). The DFN model assumes that the rock mass matrix is impermeable and the groundwater only flows in the fractures (Nejadi et al. 2017). This method has been widely used in academia because it clearly explains the structural characteristics of seepage in the fractured rock mass and is closer to the actual situation (Zhang et al. 2019; Lopes et al. 2022). However, in practical engineering, it is difficult to accurately establish the DFN model containing all fractures, and it also needs a huge calculation workload (Baghbanan & Jing 2007; Hekmatnejad et al. 2018). Therefore, this DFN model is difficult to popularize in practical engineering and large-scale hydrogeological problems (Huang et al. 2021). Most of the other models are based on these two basic models, such as the dual-porosity model (Jerbi et al. 2017), the dual-porosity and dual-permeability model (Zhang et al. 2021; Huang et al. 2022), the triple-porosity and dual-permeability model (Wang et al. 2019; Hu et al. 2020), and the multi-porosity model (Yan et al. 2016).

In analyzing many hydrogeological and engineering geological problems, the EPM model is the most practical and effective method to study the permeability characteristics of the FGM (Scanlon et al. 2003; Jaunat et al. 2016; Mohtashami et al. 2022). However, due to its simplifications of the fracture network, including overlooking the complexities of fracture distribution, connectivity, and shape, when considering the distribution of water head and flow velocity, the EPM model used in numerical calculations sometimes cannot reflect the real seepage state of the water flow field in the FGM (Zeng et al. 2021). It could result in unacceptable errors in the simulation results of the EPM model, which influence the decision-making of engineering geological problems (Khoei et al. 2016; Yan et al. 2019). Therefore, to ensure the reliability of engineering decisions, it is crucial to calculate the error of the EPM model to determine whether the simulation result of the model is acceptable, providing a scientific basis for evaluating the suitability of the EPM model more accurately.

The main influencing factors of the permeability of the FGM are the geometric characteristics of the fracture, such as fracture length, aperture, surface roughness, direction, density, and filling (Li et al. 2021). Many researchers found that natural fractures show fractal characteristics, and the geometric characteristics of fractures are related to the fractal dimension (Erol et al. 2017; Deng & Zhu. 2020). Thus, the permeability characteristics of the FGM are related to the fractal dimension of fracture (Li et al. 2019; Feng et al. 2020). How to employ the fractal dimension to predict the permeability of a fracture network has been an interesting issue (Huang et al. 2017; Liu et al. 2017). Some scholars have put forward some analytical and empirical models to predict the permeability of fracture networks and predict the representative elementary volume (REV) size of FGM permeability (Jafari & Babadagli 2012; Zhang et al. 2018). However, the influence of the fractal dimension of fracture on the suitability of the EPM model is still unclear, and it also lacks a prediction model for quantitatively evaluating the suitability of the EPM model under different fractal dimensions.

The aim of this study is to look at the effect of the fractal dimension and the length–width ratio of the simulation domain on the suitability of the EPM model. Using the Monte Carlo method, fracture network models with specified fractal dimensions and length–width ratios are generated. The EPM and DFN models are subsequently employed to simulate seepage, with the DFN model providing the reference water head values in the simulation domain. This enables a quantitative analysis of the EPM model's accuracy by comparing the water head at corresponding nodes with those predicted by the DFN model. The study concludes with an evaluation of the EPM model's suitability, exploring the influence of the fractal dimension and the length–width ratio on it. A new empirical formula is proposed to calculate the suitability of the EPM model based on these geometric parameters, providing a criterion for model selection in hydrogeological simulations. This approach introduces a novel methodology for assessing the practical suitability of the EPM model, enhancing its applicability in addressing engineering geological challenges.

After the number of fractures is updated from N_t to N'_t, the fractal length distribution of fractures in any L_n of the DFN model can be generated by formula (2). Thus, a series of stochastic DFN models with different fractal dimensions can be generated by the Monte Carlo method.

In this study, the Monte Carlo method is used to simulate a two-dimensional (2D) DFN. The Monte Carlo method is also called the random sampling technique or the statistical experimental method based on the law of large numbers in probability. At present, the Monte Carlo method is generally considered a relatively accurate method (Liu et al. 2019). The main means of the Monte Carlo method is to obtain the random numbers for statistical experiments, and the generation of random numbers that conform to a certain probability distribution function is the basis of applying the Monte Carlo method. According to the probability distribution form and numerical characteristics of various geometric parameters of fractures, the Monte Carlo method can be used to simulate the fracture network model. Based on the Monte Carlo method, a fracture network generation program is compiled in this study. This program can generate a fracture network that characterizes the fracture information of the rock mass and store the generated fracture network model data in the form of a matrix for subsequent calculation. The Monte Carlo method is a reasonable, feasible, and effective approach at the current research level of fracture investigation and simulation (Gao et al. 2016; Guo et al. 2020).

The fractures in the DFN model are generally regarded as parallel plates, in which the flow rate of the fracture surface is proportional to the cubic of the fracture width. The flow pattern of the fluid flow in fractured rock mass is generally linear and obeys the cubic law (Chen et al. 2015). Therefore, laminar flow without considering turbulence is used in this study. In tight rock masses such as granite rock and basalt, the permeability of the rock matrix is far less than that of fracture (Wang et al. 2016). Therefore, in this study, the permeability of the rock matrix is ignored, and the water flow in fractured rock mass only occurs in fractures.

In order to quantify the influence of the fractal dimension and the length–width ratio of the simulation domain on the suitability of the EPM model, the suitability of the EPM model is solved by comparing simulation results of the EPM model with those of the DFN model, and the simulation results of the DFN model are taken as the reference value of water head in the study area because of the higher accuracy of the DFN model in characterizing water head distribution. The head errors of each node are calculated by comparing the water heads of the EPM model with those of the DFN model on corresponding internal nodes. Three kinds of error are calculated: root mean square error (RMSE), mean absolute error (MAE), and relative error (RE).

The above two indicators are used to describe the error between the water head distribution of the EPM model and that of the DFN model quantificationally. The difference between them is that RMSE squares the deviation once. If the maximum deviation is large, RMSE will be amplified.

By changing the fractal dimension and the ratio of length to width of the simulation domain, this study analyzes their influence on the suitability of the EPM model and judges whether the EPM model can be used in numerical or analytical calculation for an FGM site. The fractal dimension is progressively increased from 1.3 to 1.8, and the length–width ratios of the model are divided into 1:2, 1:1.5, 1.5:1, and 2:1. Then, the Monte Carlo method is used to generate 24 groups of different fracture networks, and each group repeats the simulation 200 times to produce 200 realizations for reducing the influence of randomness on the results. Finally, 4,800 numerical simulations are conducted. The 24 groups of the fracture network with different fractal dimensions and length–width ratios of the simulation domain are simulated by the EPM model and the DFN model, respectively, and the water head of each node in the fracture network is calculated, and the simulation results of the two models are compared.

As can be seen from the results, the water head distribution results can be divided into three types. First, the simulated water head of the DFN method is evenly distributed on both sides of the line (Figures 5(f), 6(d), 6(f), 7(a), 7(b), 7(e), 8(a), 8(e), and 8(f)). In this case, the deviation between the simulation results of the EPM method and the DFN method is small, and the EPM method is more suitable. Second, the distribution of simulated water heads by the DFN method is on one side of the line (Figures 5(a), 5(c), 5(d), 6(c), 7(d), and 7(f)), which indicates that the simulation results by the EPM method are larger or smaller than those by the DFN method. At the left and right boundaries of the simulation domain, the deviation of the simulation results by the EPM method is smaller and the accuracy of the simulation results is higher, while in the middle of the simulation domain, the deviation of the simulation results is larger. Finally, the DFN simulation results show that the part of the area is above the EPM method simulation results and the other part is below the EPM method simulation results (Figures 5(b), 5(e), 6(a), 6(b), 6(e), 7(c), 8(b), 8(c), and 8(d)), which indicates that the EPM model simulation results are better near the middle position of the model, and the deviation of other positions is larger. Moreover, for this kind, the simulation results tend to show that the water head at positions close to the constant head boundary is more inclined to the boundary head value.

Through comprehensive analysis, it is found that there is a strong correlation between the distribution pattern of the water head in the simulation domain and the location distribution of fracture intersections. When the intersections are relatively evenly distributed in the simulation domain (Figures 5(f), 6(d), 6(f), 7(a), 7(b), 7(e), 8(a), 8(e, and 8(f)), the permeability difference at different positions in the simulation domain is small, and the water head results show the first kind of distribution form. When the number of intersections in the middle of the simulation domain is moderate (Figures 5(a), 5(c), 5(d), 6(c), 7(d), and 7(f)), the permeability of the middle of the simulation domain is higher, and the water head results often show the second kind of distribution form. When the number of intersections in the middle of the simulation domain is fewer than that at the left and right boundaries of the simulation domain (Figures 5(b), 5(e), 6(a), 6(b), 6(e), 7(c), 8(b), 8(c), and 8(d)), the permeability of the simulation domain near the left and right boundaries is better, but the water head results often show the third kind of distribution form. The simulation results show that there are horizontal points of water head distribution in Figures 5–8, which show that the water heads of successive positions in the study area are equal, indicating that part fractures are connected with other fractures with one side but the other side is a dead end. There are also some water head distribution results in some groups with larger slopes, which indicate that the number of fractures in this position is small and the permeability is small in the study area.

The water head distribution of the DFN model is used as the reference, and the water heads at the same position of the EPM model and the DFN model are compared. The changing trends of water head difference with fractal dimensions and different length–width ratios are similar. It can be seen that the difference in water heads between the two models is greatest when the D_f is 1.3 (Figures 5(a), 6(a), 7(a), and 8(a)). When the D_f is higher than 1.3, with the increase of the fractal dimension, the difference in water head between the EPM model and the DFN model simulation results at the same position shows an increasing trend, and the distribution of scatter points becomes more dispersed according to the indication of R². However, there are some distinctions between the variation trend of differences in water head distribution with different length–width ratios. When the length–width ratios are 1:2 and 1:1.5, the opposite trend appears and the difference in water head gets smaller when the D_f is 1.8. When the length–width ratios are 1.5:1 and 2:1, the same phenomenon appears when the D_f is 1.7 (Figures 7 and 8).

Based on the summary and analysis of 4,800 groups of numerical simulation results, two empirical formulas are proposed to calculate the suitability of the EPM model, so as to judge whether the EPM model is suitable for a study area. According to the studies, the suitability of the EPM model can be assessed in practice. First, four constant parameters are calculated according to the length–width ratio of a study area. Then, the fractal dimension effect equation is obtained according to the four constant parameters. Finally, the fractal dimension of fractures in the study area is investigated, and the final suitability of the EPM model can be calculated by the fractal dimension effect equation.

This research used a 2D steady-state flow model. The accuracy of the EPM model's simulation results was quantitatively analyzed by comparing water head values at corresponding nodes with the DFN model. The suitability of the EPM model was evaluated, exploring the impact of fractal dimensions and length–width ratios. Based on these geometric parameters, an empirical formula to calculate the suitability of the EPM model was proposed. This method inherently limits our ability to encapsulate the full spectrum of fluid dynamics within fractured geological structures. Moreover, it only includes comparisons of water head values, does not compare flow velocities, and does not adequately capture the scale-specific effects of fractures. Although this study provides valuable insights into the suitability of the EPM method under specific conditions, it overlooks key geometric parameters of fractures, such as fracture roughness and the hydraulic properties of the rock matrix. These parameters could also influence the effectiveness of the EPM method; hence, further investigation is necessary in subsequent studies to enhance the model's accuracy and extend its applicability to a broader range of geological settings.

In this study, a DFN generator based on the fractal dimension and Monte Carlo methods is compiled. Based on the DFN generator, 24 groups of models with different relevant parameters are constructed, and each group of models generates 200 realizations to reduce the influence of the randomness of the Monte Carlo method on results. A total of 4,800 numerical flow simulations are conducted. Based on the results of these numerical simulations, the influences of the fractal dimension of fracture and the length–width ratio of the simulation domain on the suitability of the EPM model are analyzed.

The results show that the suitability of the EPM model decreases at first and then increases with the increase of the fractal dimension under different conditions of the length–width ratio of simulation domains. When the length–width ratios are 1:2 and 1:1.5, the inflection points of suitability curves appear at D_f = 1.7, and when the length–width ratios are 1.5:1 and 2:1, the inflection points of suitability curves appear at D_f = 1.6. The same method is used to analyze the influence of the length–width ratio of the simulation domain. The influence of the length–width ratio on the suitability of the EPM model can be fitted with straight lines, and the slope of the straight lines gradually increases with the increase of D_f. Based on the analysis results of the fractal dimension and the length–width ratio of the simulation domain, both of them affect the suitability of the EPM model, and the length–width ratio has less influence on the suitability of the EPM model than the fractal dimension. An empirical formula is obtained based on the 4,800 numerical simulation results to quantify the suitability of the EPM model to the FGM in a study area. When using the empirical formula, firstly, the fractal dimension of fractures in a study area should be calculated according to fracture investigation. Then, four constant parameters are calculated according to the fitting results of the length–width ratio of the simulation domain, and the formula for calculating the suitability of the EPM model is formed. Finally, the suitability value can be obtained by the empirical formula of the suitability of the EPM method.

We have noticed that some important problems need to be solved in future research. First of all, our research is only restricted to a 2D fracture network, which will be developed into three-dimensional (3D) in the future. In 3D space, the influencing factors of the suitability of the EPM method to the FGM in 3D space are expected to be more complicated. The fractal dimension and length–width ratio can basically represent the spatial characteristics of fractures, which can be regarded as the main influencing factors of the EPM method suitability. However, there are some other geometric parameters of fractures that may have significant influences on the suitability of the EPM method that has not been comprehensively considered, such as fracture roughness and hydraulic properties of rock matrix, which will be further explored in the future. Furthermore, we plan to integrate our proposed evaluation method with fractal dimension measurement techniques to develop a comprehensive and convenient on-site assessment method. This approach is aimed at precisely determining the suitability of the EPM model in practical field settings, thereby enhancing its utility and effectiveness in real-world applications.

This study is supported by the National Natural Science Foundation of China (Nos U2267218 and 42072276). The authors would like to acknowledge the editors and anonymous reviewers for their invaluable discussion and suggestions.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.