Abstract
In this study, physical experiments and numerical simulations are combined to provide a detailed understanding of flow dynamics in fracture network. Hydraulic parameters such as pressure head, velocity field, Reynolds number on certain monitoring cross points, and total flux rate are examined under various clogging conditions. Applying the COMSOL Multiphysics code to solve the Navier-Stokes equation instead of Reynolds equation and using the measured data to validate the model, the fluid flow in the horizontal 2D cross-sections of the fracture network was simulated. Results show that local clogging leads to a significant reshaping of the flow velocity field and a reduction of the transport capacity of the entire system. The flow rate distribution is highly influenced by the fractures connected to the dominant flow channels, although local disturbances in velocity field are unlikely to spread over the whole network. Also, modeling results indicate that water flow in a fracture network, compared with that in a single fracture, is likely to transit into turbulence earlier under the same hydraulic gradient due to the influence of fracture intersections.
1. Introduction
Fractured flow in fractured rocks has received considerable attention for its significance in the fields such as radioactive waste disposal and subsurface contaminant transport analysis. In particular, the fracture flow behavior is a primary concern because of its dominant transmissivity compared with the rock matrix. The void space of rock fracture network is composed of individual fractures and intersections between those fractures, thus introducing heterogeneities in three dimensions across the network. An abrupt change in void geometry of the flow network will break the balance between pressure gradient, velocity, and flow exchange amount, which will consequently influence solute transport. Clogging in the fracture space can be easily found in natural fracture fields. Generally, transport in fractures is influenced by Taylorian dispersion and macro dispersion due to different solute velocity and various fracture apertures. The clogging-induced aperture decreasing will result in an exponential increase in macro dispersion [1]. Fracture clogging occurs when biological materials or translocated debris narrows or even blocks the effective aperture of the fracture. In the remediation of groundwater resources, biological clogging is usually used as a kind of biobarrier to control the polluted groundwater movement [2, 3]. However, clogging phenomenon may become a serious problem in operating underground storage caverns of hydrocarbon.
Study on the geological media clogging has made great progress in the past decades, among which most studies have focused on porous media, while only a few have explored the clogging of fractured media [4–7]. What is more, among the existing experimental studies of fracture clogging, research was mainly focused on the biological clogging mechanism in a single fracture, but the clogging effect on hydrodynamics of fracture networks is still unclear [8, 9]. Compared to the porous media or single fracture, flow in a fracture network has a much more complex pattern which presents a challenge to determining the flow characteristics. The relationship between the hydraulic conductivity () and the clogging was frequently used to present the effect of clogging process because clogging results in a reduction in [2]. As is a macroaveraged parameter that reflects the ease with which water can move through fractures, scientists are unable to get a deeper understanding of flow details in a fracture network system under mesoscale level. In a natural fracture system, flow generally occurs as a channeling effect that appears to be a preferential flow phenomenon governing flow and transport due to the heterogeneity [10–12]. The clogging of the small fractures and fissures may be very different from the clogging effect on the dominant factures; thus using value to determine the clogging effect may be inaccurate for understanding real flow phenomenon in fracture systems.
The spatial scales as well as the air-water flow environment can easily allow the formation of clogging in relatively small-aperture fractures, and, also, its apertures give a much higher potential for clogging during geological or biological activity processes, such as rainfall induced sediment deposition, curtain grouting in bedrock, and plant roots growing. In addition, smaller apertures are more difficult to be detected in the field and may be ignored during the conceptualization for a numerical model. Also, evidences have shown that tiny features in geological media can influence large-scale flow through an aquifer [13]. So far, it is still unclear about the impacts of clogging on hydrodynamics in a fracture network, which may affect the simplification of fracture models for application. The above discussion raises questions we hope to address. How will clogging fracture affect the flow hydrodynamics in a fracture network system? What are the contributions of small-aperture clogging to the entire flow field?
Here we explore the potential impact of clogging of an individual fracture on the main channel flow paths and the conveyance of the network by (1) characterizing the clogging fracture on different places via velocity field on certain sections; (2) experimentally testing the changes in total flux of the entire network and pressure heads on certain intersection points; and also (3) numerically examining the dominant flow paths and Reynolds number with various types of clogging conditions.
2. Experiment
2.1. Experimental Setup
We constructed a simple fracture network made by plexiglas for physical experiment studies. Figure 1 is a schematic diagram of the experimental setup. This apparatus consists of two subsystems: one is the physical fracture network and the measuring tools; the other one is the water circulating system for maintaining constant head both upstream and downstream. For the physical model, in order to examine the clogging influences, it was designed into two regions according to the different grid forms, named R1 and R2, respectively (Figure 2). The aperture value of wide channels of both R1 and R2 is 5.5 mm with fracture lengths 282.8 mm and 223.5 mm, respectively. The aperture value of the narrow channels of both R1 and R2 is 1.85 mm with fracture lengths 200 mm and 223.5 mm, respectively. The height of each fracture of both R1 and R2 is 100 mm (see Table 1). There are 85 fractures and 52 intersections in total in the entire network. The cross angle is the same on each individual region but different between the two regions, thus introducing large heterogeneity into the entire network.
Figure 3 shows the pictures of the experimental setup. A wooden horizontal platform with dimensions of 200 × 80 × 20 cm3 was first laid out. Each fracture was made with two vertical plexiglas boards. The aperture of fractures was maintained with extreme care during the model construction considering its importance to the accuracy of the experimental results because the flux is proportional to third power of the aperture. To solve this problem, copper sheets with known thickness were inserted along the upper side and down side of the aperture and kept in good contact with those two vertical planes. The two vertical planes representing the fracture are then fixed and sealed from outside using silica gel. Confined flow is studied here so both the upper side and down side of each fracture are sealed. All the fractures are then welded together to form the network on the intersection carefully. No leakage was observed during the experimental period.
(a)
(b)
(c)
(d)
Although there are numerous literatures related to fracture network flow and solute transport study, they are rarely combined with the physical experimental method. The reason is mainly due to the deficiency of current testing techniques. Indoor fracture networks that were reconstructed by artificial concrete blocks or natural rock debris can be much closer to the natural state of fractures, but it is hard to obtain the flow characteristics due to the opaque walls [14, 15]. Therefore, we constructed the physic model with transparent plexiglas. Different from orthogonal fracture regime reported in the past, the fractures in our physical model show different length, orientation, density, and aperture, which is apt to form obvious preferential flow phenomena.
2.2. Experimental Methods
A series of flow tests under different pressure heads was carried out using the physical model. The flow direction is from left to right parallel to the x-axis. Rotor flowmeters were used to measure the flux rate and the values were verified by the volumetric method. Piezometers were used to measure pressure heads at the center of each intersection. Since the hydraulic gradients along the fractures were very low in the physical model and fluctuations of piezometric heads were mostly lower than 2 mm, it would be inaccurate to measure the hydraulic heads using traditional methods with piezometric tube only. To solve this problem, a moveable laser scanning detector with a resolution of 0.1 mm was applied to measure the water-air interface in the piezometric tubes. In addition, to verify and adjust pressure in piezometric tubes, electronic manometers (HM28 series with precision of ±0.2%, Nanjing Helm Sci-tech Co., Ltd) with data logger were used to provide reference pressure head values. Although the synthetic fracture network is still much simpler than a natural fractured field, the physical fracture network composed of four groups of inclination and two different apertures can better approximate natural conditions than previous single-fracture models [16, 17]. The results of this study can help us gain insights of flow behaviors under clogging conditions.
Flow tests have been performed in the fracture network model with different hydraulic gradients. The following two sets of experiments have been carried out to fulfill different research objectives.(1)Seven sets of experiments with different hydraulic gradients have been conducted to calibrate the numerical model; see Table 2.(2)Six sets of physical experiments have been conducted to test the clogging effects on the preferential flow paths, fracture flow hydraulics, and total flux. Here we only selected a limited number of fractures around the main channels for clogging analysis. The main channels can be obviously observed through red dye tracer. The clogging cases were represented in Table 3. Fractures 21, 23, and 25 are in the same straight line, and they occupy the center line of R1 region. Fractures 23 and 41 bridge the two wide-aperture fracture channels of f36-f27-f18-f9 and f37-f28-f19-f10. Fractures 21 and 25 are located upstream and downstream of fracture 23, respectively. Fractures 46 and 48 are the two links between regions R1 and R2.
The monitoring items for each clogging case were described in Table 4. For the monitoring fractures, fracture 18 is on the main channel of R1 while fractures 65 and 79 are on the main channel of R2. Fractures 19 and 28 are in the same straight line and intersect fracture 23 at cross point 16. Fracture 76 is far from the main channel compared with other monitoring fractures. For the monitoring cross points, CP18 is the center point of the network, CP16 and CP21 are on the main channel of R1, while CP45 is on the main channel of R2, and CP31 and CP32 are the link points of R1 and R2.
We treat the clogged fractures as a kind of fully choking condition instead of a gradually blocking process that would occur in natural world. Evidences have shown that bioactivities can result in a two-order magnitude decrease in hydraulic conductivity [2]. So it makes sense to block a fracture completely in our experiments.
3. Modeling
Due to the difficulty of measuring detailed fracture flow features in the physical model, a numerical modeling approach was adopted for further analysis. The COMSOL Multiphysics model was applied in this study to simulate flow dynamics in the fracture network.
3.1. Numerical Methods
(1) Governing Equations. The Reynolds equation (i.e., local cubic law) which is simplified from the Navier-Stokes equations is commonly used to describe the fluid flow in rock fractures when the Reynolds number is low and the fracture geometry does not vary too abruptly. The Reynolds equation is expressed as where the term is usually called the fracture transmissivity and is the aperture of the idealized parallel smooth fracture.
In the past decades, the validity of the Reynolds equation and local cubic law for rock fractures have been investigated widely by using artificial or natural rock fractures [18, 19]. The common understanding is that their applicability can be guaranteed only when the Reynolds number is small (where viscous forces dominate the inertial forces) and aperture does not change abruptly [20, 21]. The deviation of the flow velocity fields from ideal parabolic profiles across the fracture aperture, which could happen when fractures are not planar and smooth as required by Reynolds equation, could have a significant impact on particle transport in rock fractures [19]. Existing literatures have concluded that the Reynolds equation overestimated the flow rate by as much as 100% and might not be suitable for estimating the flow in rock fractures [18]. To better quantify the flow properties in fracture networks, the NS equations combined with the continuity equation need to be solved directly: where is the fluid density (kg m−3), is the dynamic viscosity (kg m−1 s−1), is velocity vector, and is fluid pressure (Pa).
(2) Reynolds Number. Flow in a single fracture is closely related to its Reynolds number, which reflects the relative importance of inertial versus viscous force. The Reynolds number, Re, for flow in a circular pipe is defined as where is the flow velocity (m s−1), is the diameter of circular pipe (m), and is the kinematic viscosity (m2 s−1). For flow in a cross-section other than the circular shape, is substituted by the hydraulic radius (). In a confined flow condition with a rectangular shape of cross-section, is where is the fracture aperture (m) and is the water depth (m), equal to the height of fracture here. Thus, the Reynolds number for flow in fractures becomes
As is much smaller than , (6) can be simplified into the following equation:
To solve NS equations (2) together with continuity equation (3), we used the commercial FEM software COMSOL Multiphysics. The numerical simulation was conducted using a steady state 2D horizontal cross-section geometry model of the fracture network. Steady state simulation was used to pursue the solution under equilibrium conditions. Triangular mesh type was implemented for the finite element flow simulation in COMSOL Multiphysics (Figure 2). The domain consisted of more than 120 000 triangular elements with node spacing of less than 0.4 mm. Simulations were tested for mesh dependence.
3.2. Model Validation
In our experiments, details of flow features such as velocity distributions within a cross-section and the vertical velocity along the fractures cannot be captured directly from the physical measurements. In view of this, modeling study is needed to further understand fracture network flow behaviors.
Figure 4 shows the plot of measured pressure head versus computed values . Values of 36 measured points across the whole network were applied for the analysis. The correlation coefficient is 0.985 and the root mean square error (RMSE) is 2.042 mm. This result confirms the high reliability of the simulation model. In addition, in order to better reflect the accuracy of numerical model, we examined the flow rate under various pressure heads and compared with the simulated values. All the testing cases and simulations were under steady state condition. In each case, the flow rate was measured no less than 5 times and the average value was used to represent the steady state flow rate. The 2D simulation resulted in fairly accurate total flow rate values from the numerical models. Due to the high ratio between the height and the aperture of each individual fracture, the vertical velocity profile approximated a uniform distribution, so that the 3D effect can be neglected, which assured the success of the 2D modeling approach. Thus, the flow rate can be calculated by with the flow rate, the fracture height, and the fracture aperture.
Both the measured flow rate and computed flow rate versus averaged hydraulic gradient were plotted in Figure 5. The calculated flow rate and measured values are very close, which verifies the effectiveness of the assumption of 2D flow. Small errors still exist. When the hydraulic gradient is less than 0.05, the calculated flow rates are slightly less than the corresponding measured values. When the hydraulic gradient is greater than 0.062, the calculated value is slightly larger than the measured value, and, with the hydraulic gradient increasing, the gap also increases. Nevertheless, the relative errors of the results obtained by numerical method related to those obtained by measurement are under 1%.
In the test, there is a certain main flow channel that is composed by several fractures in the fracture network. This main channel occupies a deterministic position in the network and accounts for the vast majority of the total flow of the entire network. The main channel, which shows the least resistance, can convey the water flux effectively from region R1 to R2. In Figure 5, with a smaller hydraulic gradient , the shows a linear relationship (the dash-dot line in Figure 5). However, with the hydraulic gradient increases, the flow rate in the fracture network is increased, and the substantially deviates from the linear relationship. The trend can be fitted by a power law function with a high goodness of fit.
4. Results and Analysis
4.1. Impact of Clogging of Small-Aperture Fractures on the Network Flow Field
In view of the fact that relatively small-aperture fractures have higher potential for being clogged under biological or geological activities, the clogging of smaller aperture fractures is examined to make a quantitative analysis. The following results show that small-aperture fractures have nonnegligible influences on the dominant channel flow, and the location of the small-aperture factures plays an important role. The clogging of small fractures can change not only the pressure head and velocity field of the network, but also the Reynolds numbers in local regions that is generally used as an indicator to distinguish (laminar or turbulent) flow regime.
(1) Velocity. Figure 6 shows the velocity distribution of six fixed cross-sections under 6 kinds of clogging cases, and each monitoring section is located at the middle point of the fracture. Also, the horizontal velocity profiles under nonclogging conditions were plotted for comparing analysis. It is clear that clogging at different locations introduces different magnitude of changes in velocity along the main channels, especially for fractures 18, 19, and 28. For example, Figure 6(a) shows that clogging of fractures 23 and 41 has significant impacts on the velocity of fracture 18; the velocity increased by 34.5% and 20.7%, respectively, compared to the values of the nonclogging case. Similarly, the clogging of both fractures numbers 23 and 41 has a significant impact on the flow rate of fracture 19. Compared to the flow rate without clogging, the flow rate increased by 23.7% and 85.2%, respectively (Figure 6(b)); and, for fracture 28 which is located in the same straight line with fracture 19, the flow rates were increased 1.67 times and 1.76 times, respectively (Figure 6(c)), by clogging fractures 21 and 41. It is interesting to note that the clogging of fracture 46 brings opposite effects on fractures 19 and 28 although they are located along the same straight line: the maximum flow rate increased by 54.8% in fracture 19 but reduced by 82.5% in fracture 28. It implies that local clogging of small fracture may introduce significant changes to the local velocity field and may greatly alter the flow rate in the main flow channels.
(a)
(b)
(c)
(d)
(e)
(f)
Contribution of small-aperture fractures to the flow velocity field depends on the relative position between the main channels and small-aperture fractures. For example, the function of fractures 21 and 25 in R1 region is completely different from fracture 23. Also, according to Figures 6(a)–6(f), the clogging of fracture 48 causes relatively big changes only to fracture 19 but little influence on other monitoring sections in either R1 or R2 regions. This indicates that local clogging will affect the local flow network only and this influence is less likely to spread to the entire flow network due to the redistribution of the flux in the well-connected fracture network.
(2) Pressure Head. Figure 7 shows the pressure head changes on monitoring cross points of the network under different clogging conditions. Besides the listed monitoring cross points in Table 4 (dark bars), pressure heads of an additional group of cross points were monitored (gray color bars). At this point, the pressure heads on monitored points substantially represent the head pressure distribution of the whole network. Except for the clogging case of fracture 21, clogging of other fractures has more or less impacts on the network pressure heads. A general rule is that the clogging of the fractures will lead to increase in water pressure heads upstream and decrease in water pressure heads downstream. This phenomenon has been reflected in Figures 7(b), 7(c), 7(e), and 7(f). In our experiments, the blockage of a single fracture leads to the pressure changes in the range of −13% to 6% compared with a nonclogging case. Compared to other cases, the pressure heads on monitored points are disturbed intensely under fracture 41 clogging condition (see Figure 7(d)). Except for the unchanged pressure heads on CP8 and CP47, pressure heads on all other points are dropped, especially for point CP21, with an obvious decrease of nearly 13%. This phenomenon can be attributed to the fact that fracture 41 is acting as a key dominant flow channel within the whole network. It implies that grasping the position of a key flow channel during a geological exploration will be effective for determining subsurface transport and will greatly reduce the amount of geological curtain works.
(a)
(b)
(c)
(d)
(e)
(f)
(3) Reynolds Numbers and Flow Regime. Figure 8 illustrates fluctuations of different monitoring points along the main channels under both clogging and nonclogging conditions. This is similar to the impacts on velocity field where the local clogging of small fractures will lead to significant changes to local . The clogging of fracture 41 made a sudden fall of the at cross point 21 (falling from 568.12 to 108.89). It implies that fracture 41 consists of a part of the main flow channels. Thus, the function of small-aperture fractures is determined by the spatial location and also the relative connectivity with main channels, but not the geometric property only.
(a)
(b)
Reynolds number is an important parameter usually used for characterizing flow regime. Although it was difficult for us to measure the critical Reynolds number for turbulence transition with a limited number of test cases, the vortex indicating flow approaching transition (from laminar flow to turbulent flow) was clearly obtained through numerical particle tracing. Figure 9 shows the particle trajectories around four typical crossing points. These particle trajectories are calculated by using Lagrangian tracking method. All the particles are treated as massless particles that follow the velocity field. Theoretically, if particles travel by strictly following the streamlines, they will not go into these trapping zones with closed streamlines that will not occur in reality. If particles jump between streamlines due to other physical processes such as molecular diffusion, such trapping might occur to reduce the particle travel speed and generate long tails in the breakthrough curves. Circles and trapping zones are often observed around the fracture intersections which imply that local rotational flow exists in these areas. The vortex is much more obvious in wider fractures compared to the flow in narrow fractures. Vortexes were often observed even when the pressure head difference was small between the upstream and downstream directions, especially in or around the main fractures.
(a) CP15
(b) CP18
(c) CP6
(d) CP37
It is unclear if it is adequate to use the observed vortexes as the criterion to justify laminar or turbulent flow in fractures; however, the observed vortexes can be considered as a foregleam of the coming turbulence. Results of the experiments showed that vortexes would occur in a fracture network under a relatively low Reynolds number (less than 900). This is different from previous results showing that fracture flow remains laminar when the Reynolds number is below 2320 [22, 23]. Similar research results have been reported in other fracture flow studies; for example, Qian et al. [16] found that the possible Reynolds number for turbulent transition in a fracture is much lower than flow in a pipe where the Reynolds number must be larger than 100,000 to make flow turbulent [24]. Although Qian et al. [16] attributed the lower Reynolds numbers to the much smaller hydraulic radius () of a fracture as compared to hydraulic radius of flow in a pipe (the diameter), it is not enough to explain the case in fracture network. The deflection flow phenomena in a fracture network make the hydraulic property totally different from flow in a single fracture. In order to illustrate the influence of fracture intersections on the hydrodynamics within a network, one main flow channel in region R1 composed of f36, f27, f18, and f9 was extracted from the network for further modeling analysis. The pressure head and velocity along the center line of fracture f36-f27-f18-f9 were compared with a single fracture without any intersections; see Figure 10. The pressure head difference between the upstream and downstream directions of both the two models is set the same. It is clear that, with intersections, the pressure head is no longer linearly distributed but shows a nonlinear pattern. Under steady flow state, the flow velocity along the fracture is no longer a constant value but shows an apparently discontinued, step-function distribution. And accordingly, the Reynolds number of each fracture within the networks is not a constant value as in a single fracture. This indicates that, even under the same hydraulic gradient, the turbulence phenomena will occur earlier in the fracture networks. With the increasing of hydraulic gradient, flow in parts of the fractures will transit into turbulent first. This finding is important for evaluating the transport properties in fractured media as the main channel flow would be highly affected by small but well-connected fissures around the channels. Ignoring those small fractures in the conceptual models may lead to biased results in prediction and evaluation of contaminant transport in fractured aquifers.
(a)
(b)
4.2. Local Clogging on the Preferential Flow Paths
In nature geological mediums, faults, joints, and micro fissures have crucial influences on hydrological properties of most geological formations. Their impact on water flow and contaminant migration has been investigated extensively because they serve as preferential flow paths [10]. Variable fracture apertures have been shown to account for water channeling and for spatially variable flow patterns within the fracture void. The concentrated flow along these nested preferential flow paths also enhances solute and contaminant transport. Here we examine the impacts of local clogging on the preferential flow paths. Figure 11 shows the flux distribution in the network under clogging and nonclogging conditions. Based on the result, the flux rate values in fractures are grouped into four levels that are represented by the width of the solid lines. The main channels are identified by the flux rate values. Preferential flow is evident in the entire network under both clogging and nonclogging conditions. In R1 region, the dominant channels are those small-aperture fractures with horizontal direction while in region R2 they mainly exist in wide fractures. Thus, for the fracture network, the flux distribution is closely related to the spatial location and orientation of each fracture, but not only related to the aperture size. This result is consistent with the conclusions by Dahan et al. [4]. In addition, comparing Figures 11(b)–11(g) with Figure 11(a), the local fracture clogging resulted in flux redistribution around the clogged fracture only but did not impose significant impact on flux rate distribution in the entire network, which is consistent with the above conclusions on the velocity field analysis.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
Figure 12 shows that the local clogging can reduce the total flux of the fracture network. The reduction of flux depends on the relative contribution of each fracture. Result of Case 1 shows that clogging of fracture 21 leads to nondiscrimination on the flux compared to nonclogging condition. This is because fracture 21 is not a portion of a main channel in the network. On the other hand, fracture 46 plays a significant role in forming the main channel and it also acts as the main link for regions R1 and R2, so that the clogging of fracture 46 resulted in significant changes of flow in the main channel, a 6.8% reduction of flux. In a certain fracture network, the preferential flow path means minimum energy loss along such a path. Clogging of a fracture along preferential flow paths will break the energy balance in a fracture system.
5. Conclusions
In summary, the following conclusions can be drawn from this study.(1)Clogging of small-aperture fracture can lead to significant changes in pressure heads and flow rate distribution in fractures adjacent to the clogged fracture. This effect is related to space position of the clogged fracture. Although in this study the effect is mainly demonstrated in cases where the clogged small fractures are well connected to the main flow channels, it indicates that local clogging of fissures around the main flow channels cannot be ignored in a conceptual numerical model.(2)Based on the small-aperture clogging test, clogging of small fissures that well connected to main flow channels only produces minor impacts on preferential flow paths. The change can be reflected by the flux at local scales and is less likely to propagate through the whole network. However, even for a well-connected fracture network, the clogging of small fractures may reduce the transmissivity of the entire network.(3)Due to the interference of the intersections in the fracture network as well as the water deviation effect, even under a steady flow condition, compared to the flow properties in a single fracture, the flow velocity along the fracture is no longer a constant value but shows a step-function pattern. Also, the pressure head distribution in the fractures no longer changes linearly but shows nonlinear properties. Finally, numerical modeling results suggest that the turbulence phenomena will occur earlier in the fracture networks than in a single fracture.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
This research was supported by the National Natural Science Foundation of China (no. 51349015 and no. 51297045).