By Gowri Srinivasan
Fractures are the primary pathways for fluid flow in otherwise low-permeability subsurface media, such as shale or granite. Applications where flow and transport through fractured media are central to informed decision-making via predictive capability have recently inspired a great deal of interest. Due primarily to the increased availability of natural gas, which produces 50 percent less carbon dioxide than coal, U.S. emissions have dropped to their lowest levels in 20 years. As a result, the production of hydrocarbons from shale formations—which involves hydrofracturing rock to establish fracture connectivity and extracting the natural gas flowing out of the fractures—is a topic of ongoing debate and research. Moreover, as countries like North Korea continue to improve their ability to conduct low-yield nuclear tests, chemical signature detection (e.g., xenon migration through fractured rock) provides the definitive smoking gun when used with conventional seismic methods.
Ignoring the topological and geometrical properties of these fracture networks and modeling subsurface fractures as a continuum with certain effective properties falls short of capturing key behaviors , including flow channeling and arrival times of particles advecting with the flow. It is thus desirable to retain the underlying structure of these fractured systems, resulting in the emergence of an alternative discrete fracture network (DFN) modeling methodology. Recognition that fracture geometry and network topology play a critical role in determining quantities of interest relating to flow and transport in fractured subsurface media is a key distinguishing factor of DFN modeling from standard continuum models.
Until recently, DFN models were limited to one-dimensional pipe-network approximations , two-dimensional systems, or relatively small three-dimensional (3D) systems . However, recent advances in high-performance computing have enabled flow and transport simulations in large, explicit 3D DFN representations . The 3D DFN software dfnWorks , which recently won an R&D 100 Award, assigns each network fracture a shape, location, aperture, and orientation by sampling distributions whose parameters are based on site characterization. Figure 1a shows a typical network consisting of 7,200 fractures of varying size, orientation, and permeability. Once meshed, this system has nearly 16 million nodes and 32 million triangular elements. Figure 1b depicts a close-up of the mesh representation; the flow equations are solved explicitly on this detailed representation. The increase in model fidelity comes at a huge computational cost because of the large number of mesh elements required to represent thousands of fractures with sizes that range over several orders of magnitude (from mm to km) and accurately resolve large pressure gradients at fracture intersections. Given that most sites can only be characterized statistically, computing uncertainty bounds for the purpose of decision-making requires thousands of runs.
One can think of a graphical representation of a fracture network, where nodes and edges inherit geophysical and geometric properties of that network, as a coarse-scalerepresentation that preserves the key topological properties . Optimal assignment of edge weights, a topic of research itself, can enhance the effectiveness of the reduced graph representation to emulate the DFN. Furthermore, the weighted graph framework seamlessly lends itself to numerically solving flow and transport equations on the nodes, resulting in efficient and accurate DFN emulators to use instead of the computationally expensive DFN model. The reduced computational burden is particularly valuable in the context of uncertainty quantification.
The mapping in Figure 2b represents each DFN fracture as a graph node and each fracture intersection as an edge connecting the nodes. This can be viewed as a projection of the bipartite representation on the fracture nodes. One can easily extract measures of network connectivity (degree) and importance of individual fractures (centrality) under this mapping. Predictions of first arrivals are critical for applications such as nuclear nonproliferation, since trace gases like xenon migrate in high concentrations early on. In this case, one can use a mapping based solely on the topology shown in Figure 2b to determine the shortest paths between source and sink through a fracture network . Figure 2c depicts a mapping that represents fractures as collections of edges and each intersection as a graph node; this is a projection of the bipartite representation onto the intersection nodes. Such mapping allows graph edges to inherit hydrological properties, such as permeability and in-plane geometry, and can simulate flow and transport in a computationally efficient manner. For DFN transport modeling, the breakthrough curve is a typical quantity of interest (QOI) that describes the mass breakthrough of the transported species as a function of time as it crosses a control plane at the outlet. The appropriate mapping depends on the QOI.
Figure 3 illustrates both the probability distribution function of arrival times at the outlet and the cumulative distribution function for a few of the relevant graph-based reduced-order models. A comparison to the high-fidelity DFN model is also shown to gauge the effectiveness of these models in emulating QOIs. The aggressive shortest path, which represents 10 percent of the entire network, can serve as a surrogate at very early times but deviates significantly beyond one percent of the total simulation time. This further underscores large-scale structure’s dominance among the factors controlling transport within the network, especially at early times. Furthermore, one can use the structure to identify flow channels. Machine learning (ML) approaches have shown initial success in reducing the network size to 25 percent by classifying fractures as participating in the flow . The pruned network predicted by ML techniques is a less aggressive reduction in fractures and tracks the DFN model well in regard to peak arrival.
In field applications such as hydraulic fracturing or environmental remediation, the entire breakthrough curve is typically of interest to measure production or remediation efficiency, respectively. The graph-based flow and transport models can run with up to five orders of magnitude of computational savings, even though the tradeoff is a systematic—but remediable—discrepancy from the mesh-based DFN solution. The corrected graph transport solution provides the closest match to the DFN solution at both early and peak arrival times, although discrepancies exist in the tail behavior. As mentioned previously, each of these methods is an acceptable surrogate for the full DFN model, and choosing the appropriate model depends heavily on the questions being asked.
It is important to note that one can always compare the aforementioned graph-based reduced models with the original higher fidelity DFN representation. The recent explosion of ML in geosciences and other fields raises the question of whether the developed emulators can predict QOIs outside of the training data, and what might be considered an extrapolation beyond the training regime. This criticism is particularly well-founded when learning from limited experimental data or already approximated models. However, 3D DFN models have been extensively used and validated against field observations, e.g., at the Äspö Hard Rock Laboratory in southern Sweden. This affords reasonable confidence that the “ground truth” furnished by the high-fidelity DFN models provides a solid foundation for additional insight into fractured systems.
Acknowledgments: The author thanks the following collaborators at Los Alamos National Laboratory for contributing to the research presented in this article: Aric Hagberg, Jeffrey Hyman, Satish Karra, Daniel O’Malley, David Osthus, Jamal Mohd-Yusof, Shriram Srinivasan, and Hari Viswanathan. The author also thanks Los Alamos National Laboratory’s Laboratory Directed Research and Development Program for support throughout this project.
 Bogdanov, I., Mourzenko, V., Thovert, J.-F., & Adler, P. (2007). Effective permeability of fractured porous media with power-law distribution of fracture sizes. Phys. Rev. E, 76(3), 036309.
 Cacas, M.C. (1990). Modeling fracture flow with a stochastic discrete fracture network: Calibration and validation: 1. the flow model. Wat. Resc. Res., 26(3), 479-489.
 Djidjev, H., O’Malley, D., Viswanathan, H., Hyman, J.D., Karra, S., & Srinivasan, G. (2017). Learning on graphs for predictions of fracture propagation, flow and transport. In 2017 IEEE International Parallel and Distributed Processing Symposium Workshops (pp. 1532-1539). Lake Buena Vista, FL.
 Hyman, J. D., Hagberg, A., Srinivasan, G., Mohd-Yusof, J., & Viswanathan, H. (2017). Predictions of first passage times in sparse discrete fracture networks using graph-based reductions. Phys. Rev. E., 96(1), 013304.
 Hyman, J.D., Karra, S., Makedonska, N., Gable, C.W., Painter, S.L., & Viswanathan, H.S. (2015). dfnWorks: A discrete fracture network framework for modeling subsurface flow and transport. Comp. Geosci., 84, 10-19.
 Painter, S., & Cvetkovic, V. (2005). Upscaling discrete fracture network simulations: An alternative to continuum transport models. Wat. Resc. Res., 41(2).
 Valera, M., Guo, Z., Kelly, P., Matz, S., Cantu, A., Percus, A.,…,Viswanathan, H. (2017). Machine learning for graph-based representations of three-dimensional discrete fracture networks. Comput. Geosci. In press.
Gowri Srinivasan is an applied mathematician and scientist in the Theoretical Division of Los Alamos National Laboratory. Her research interests include reduced-order models, uncertainty quantification, and machine learning.