Machine learning for graph-based representations of three-dimensional discrete fracture networks

Abstract

Structural and topological information play a key role in modeling flow and transport through fractured rock in the subsurface. Discrete fracture network (DFN) computational suites such as dfnWorks (Hyman et al. Comput. Geosci. 84, 10–19 2015) are designed to simulate flow and transport in such porous media. Flow and transport calculations reveal that a small backbone of fractures exists, where most flow and transport occurs. Restricting the flowing fracture network to this backbone provides a significant reduction in the network’s effective size. However, the particle-tracking simulations needed to determine this reduction are computationally intensive. Such methods may be impractical for large systems or for robust uncertainty quantification of fracture networks, where thousands of forward simulations are needed to bound system behavior. In this paper, we develop an alternative network reduction approach to characterizing transport in DFNs, by combining graph theoretical and machine learning methods. We consider a graph representation where nodes signify fractures and edges denote their intersections. Using random forest and support vector machines, we rapidly identify a subnetwork that captures the flow patterns of the full DFN, based primarily on node centrality features in the graph. Our supervised learning techniques train on particle-tracking backbone paths found by dfnWorks, but run in negligible time compared to those simulations. We find that our predictions can reduce the network to approximately 20% of its original size, while still generating breakthrough curves consistent with those of the original network.

Appendices

Appendix: Flow equations and transport simulations

We assume that the matrix surrounding the fractures is impermeable and there is no interaction between flow within the fractures and the solid matrix. Within each fracture, flow is modeled using the Darcy flow. The aperture within each fracture is uniform and isotropic, but they do vary between fractures and are positively correlated to the fracture size [21]. Locally, we adopt the cubic law [52] to relate the permeability of each fracture to its aperture. We drive flow through the domain by applying a pressure difference of 1MPa across the domain aligned with the x axis. No flow boundary conditions are applied along lateral boundaries and gravity is not included in these simulations. These boundary conditions along with mass conservation and Darcy’s law are used to form an elliptic partial differential equation for steady-state distribution of pressure within each network

$$ \nabla \cdot (k \nabla P ) = 0~. $$

(1)

Once the distribution of pressure and volumetric flow rates are determined by numerically integrating (1) the methodology of Makedonska et al. [33] and Painter et al. [41] and are used to determine the Eulerian velocity field u(x) at every node in the conforming Delaunay triangulation throughout each network.

The spreading of a nonreactive conservative solute transported is represented by a cloud of passive tracer particles, i.e., using a Lagrangian approach. The imposed pressure gradient is aligned with the x axis and thus the primary direction of flow is in the x direction. Particles are released from locations in the inlet plane x₀ at time t = 0 and are followed until the exit the domain at the outlet plane x _L The trajectory x(t; a) of a particle starting at a on x₀ is given by the advection equation, $\dot {\mathbf {x}(t; \mathbf {a})} = \mathbf {v}(t;\mathbf {a})$ with x(0; a) = a where the Lagrangian velocity v(t; a) is given in terms of the Eulerian velocity u(a) as v(t; a) = u[x(t; a)]. The mass represented by each particle m(a) and the breakthrough time at the outlet plane, τ(x _L ; a) of a particle that has crossed the outlet plane, x _L = (L,y,z) is can be combined to compute the total solute mass flux ψ(t) that has broken through at a time t,

$$ {\Psi}(t,\mathbf{x}_{L} ) = \frac{1}{M} \int\limits_{{\Omega}_{a}} d m(\mathbf{a}) \delta[t - \tau(\mathbf{x}_{L},\mathbf{a})]~, $$

(2)

where Ω _a is the set of all particles. Here mass is distributed uniformly amongst particles, i.e., resident injection is adopted. For more details about the injection mode see Hyman et al. [25].

Algorithm description

2.1 Random forest

The random forest method is based on constructing a collection of decision trees. A decision tree [44] is a tree whose interior nodes represent binary tests on a feature and whose leaves represent classifications. An effective way of constructing such a tree from training data is to measure how different tests, also called splits, separate the data. The information gain measure compares the entropy of the parent node to the weighted average entropy of the child nodes for each split. The splits with the greatest information gain are executed, and the procedure is repeated recursively for each child node until no more information is gained, or there are no more possible splits. A limitation of decision trees is that the topology is completely dependent on the training set. Variations in the training data can produce substantially different trees.

The random forest method [17, 18] addresses this problem by constructing a collection of trees using subsamples of the training data. These subsamples are generated with replacement (bootstrapping), so that some data points are sampled more than once and some not at all. The sampled “in-bag” data points are used to generate a decision tree. The “out-of-bag” observations (the ones not sampled) are then run through the tree to estimate its quality [30]. This procedure is repeated to generate a large number (hundreds or thousands) of random trees.

To classify a test data point, each tree “votes” for a result. This provides not only a predicted classification, determined by majority rule, but also a measure of certainty, determined by the fraction of votes in favor. The use of bootstrapping effectively augments the data, allowing random forest to perform well using fewer features than other methods. The category with more votes is assigned to the new observation. The idea of random decision forests originated with T. Ho in 1995. Ho found that forests of trees partitioned with hyperplanes can have increased accuracy under certain conditions [17]. In a later work [18], Ho determined that other splitting methods, under appropriate constraints, yielded similar results.

Additionally, random forest provides an estimate of how important each individual feature is for the class assignment. This is calculated by permuting the feature’s values, generating new trees, and measuring the “out-of-bag” classification errors on the new trees. If the feature is important for classification, these permutations will generate many errors. If the feature is not important, they will hardly affect the performance of the trees.

2.2 Support vector machines

Support vector machines (SVM) use a maximal margin classifier to perform binary classification. Given training data described by p features, the method identifies boundary limits for each class in the p-dimensional feature space. These boundary limits, which are (p − 1)-dimensional hyperplanes, are known as local classifiers, and the distance between the local classifiers is called the margin. SVM attempts to maximize this margin, making the data as separable as possible, and defines the classifier as a hyperplane in the middle that separates the data into two groups. The data points on the boundaries are called support vectors, since they “support” the limits and define the shape of the maximal margin classifier.

Formally, a support vector machine attempts to construct a hyperplane,

$$ \langle \mathbf{a}, \mathbf{x}_{i} \rangle + b = 0~, $$

(3)

that partitions data points x _i into disjoint sets, for i ∈{1,…,n}. To define the hyperplane we must determine the unknown coefficients a = (a₁,…,a _p ) and b. SVM seeks to determine a hyperplane that separates the two sets A and B and leaves the largest margin.

Let 〈a,x _i 〉 + b = ± 1 be the normalized equations for the boundaries of the hyperplane margin. Data points x _i on either side of the margins of the hyperplane will lie in either A and B depending upon whether they satisfy either

$$ \langle \mathbf{a}, \mathbf{x}_{i} \rangle + b > 1\qquad {\text or} \qquad \langle \mathbf{a}, \mathbf{x}_{i} \rangle + b < -1~. $$

(4)

Define an indicator function S(x) by

$$ S(\mathbf{x}) =\left\{\begin{array}{l} 1 \qquad \mathbf{x} \in A\\ -1 \qquad \mathbf{x} \in B \end{array}\right. $$

(5)

and set S _i = S(x _i ). Then we can combine (4) into a single inequality (〈a,x _i 〉 + b)S _i ≥ 1 for all x _i ; with equality holding for support vectors, the nearest points to the margin.

In most cases the sets A and B may only be close to linearly separable. To account for this possibility we introduce slack variables ξ _i ≥ 0,

$$ (\langle \mathbf{a}, \mathbf{x}_{i} \rangle + b) S_{i} \ge 1 - \xi_{i} \qquad \forall~i~, $$

(6)

that allows x _i corresponding to ξ _i > 0 to be incorrectly classified with the ξ _i being used as a penalty term. Distances ρ₁ and ρ₂ from the coordinate origin to the margin boundary are given by ρ₁ = −(b + 1)/∥a∥ and ρ₂ = −(b − 1)/∥a∥ where ∥a∥ is the Euclidean length of a. The margin width is the distance between these two lines d = ρ₂ − ρ₁ = 2/∥a∥. We seek to maximize d or, equivalently, minimize ∥a∥, over the set of training data (y₁,…,y _m ), subject to the linear constraints (6). Typically one works with the Lagrangian formulation of the constrained optimization problem. The dual Lagrangian problem is to minimize the objective function

$$\begin{array}{@{}rcl@{}} L(\mathbf{a}, b, \boldsymbol{ \xi}, \gamma, \boldsymbol{\delta}) &\,=\,& \frac{1}{2}{\| \mathbf{a}\|}^{2} \,-\, {\sum}_{j = 1}^{m} \gamma_{j} [ (\langle \mathbf{a}, \mathbf{y}_{j} \rangle \,+\, b) S_{j} \,-\, 1 + \xi_{j}] \\ && + C {\sum}_{j = 1}^{m} \xi_{j} - {\sum}_{j = 1}^{m} \delta_{j} \xi_{j}~ \end{array} $$

(7)

subject to the non-negativity constraints of the Lagrange multipliers γ _j ,δ _j , and ξ _j ≥ 0 and obtain a and b. The penalty, or margin parameter, C is a regularization term that controls how many points are allowed to be mislabeled in the SVM hyperplane construction; smaller values of C allow for more points to be mislabeled. A solution to this optimization problem defines a and b. The SVM classifier is then given by the sign of the decision function,

$$ \chi(\mathbf{x}_{i}) = \text{sgn} (\langle \mathbf{a}, \mathbf{x}_{i} \rangle + b). $$

(8)

SVM falls into the category of kernel methods, a theoretically powerful and computationally efficient means of generalizing linear classifiers to nonlinear ones. For instance, on a two-dimensional surface (p = 2 features), instead of the line described by Eq. 3, we can choose a polynomial curve or a radial loop. Equation 8 may then be written in the form

$$ \chi(\mathbf{x}_{i}) = \text{sgn}\left( \alpha_{0} + {\sum}_{j = 1}^{m}\alpha_{j} K(\mathbf{x}_{i},\mathbf{y}_{j})\right), $$

(9)

where K(x _i ,y _j ) is an appropriate kernel function of x _i and training point y _j . Radial kernels often provide the best classification performance, but at higher computational costs.