A Computational Workflow for Flow and Transport in Fractured Porous Media Based on a Hierarchical Nonlinear Discrete Fracture Modeling Approach

: Modeling ﬂow and transport in fractured porous media has been a topic of intensive research for a number of energy- and environment-related industries. The presence of multiscale fractures makes it an extremely challenging task to resolve accurately and e ﬃ ciently the ﬂow dynamics at both the local and global scales. To tackle this challenge, we developed a computational workﬂow that adopts a two-level hierarchical strategy based on fracture length partitioning. This was achieved by specifying a partition length to split the discrete fracture network (DFN) into small-scale fractures and large-scale fractures. Flow-based numerical upscaling was then employed to homogenize the small-scale fractures and the porous matrix into an equivalent / e ﬀ ective single medium, whereas the large-scale fractures were modeled explicitly. As the e ﬀ ective medium properties can be fully tensorial, the developed hierarchical framework constructed the discrete systems for the explicit fracture–matrix sub-domains using the nonlinear two-point ﬂux approximation (NTPFA) scheme. This led to a signiﬁcant reduction of grid orientation e ﬀ ects, thus developing a robust, applicable, and ﬁeld-relevant framework. To assess the e ﬃ cacy of the proposed hierarchical workﬂow, several numerical simulations were carried out to systematically analyze the e ﬀ ects of the homogenized explicit cuto ﬀ length scale, as well as the fracture length and orientation distributions. The e ﬀ ect of di ﬀ erent boundary conditions, namely, the constant pressure drop boundary condition and the linear pressure boundary condition, for the numerical upscaling on the accuracy of the workﬂow was investigated. The results show that when the partition length is much larger than the characteristic length of the grid block, and when the DFN has a predominant orientation that is often the case in practical simulations, the workﬂow employing linear pressure boundary conditions for numerical upscaling give closer results to the full-model reference solutions. Our ﬁndings shed new light on the development of meaningful computational frameworks for highly fractured, heterogeneous geological media where fractures are present at multiple scales.


Introduction
Modeling fluid flow and transport in fractured porous media on the field scale, both accurately and efficiently, has been a long-standing issue in a number of energy-and environment-related disciplines, such as production of oil and gas from naturally fractured reservoirs, extraction of geothermal energies, disposal of nuclear waste underground, sequestration of carbon dioxide in deep saline reservoirs, conditions that gives rise to diagonal equivalent/effective permeability tensors, which yields some convenience for two-point flux approximation (TPFA)-based simulators. In practice, a full permeability tensor with non-zero off-diagonal elements may be more appropriate to represent the flow properties of the rock matrix together with the small-scale fractures as the DFNs often have a predominant orientation because of the state of the stress field when the fractures were formed.
In this work, we follow similar ideas presented in [16] and build a computational workflow for modeling flow and transport in fractured porous media. The aim of this research was to investigate and compare the effect of different boundary conditions on the upscaling of small-scale fractures. Specifically, we considered the constant pressure drop boundary condition that gives rise to a diagonal equivalent/effective tensor and the linear pressure boundary condition that leads to a full tensor with non-zero off-diagonal elements. Once the small-scale fractures were upscaled, the large-scale fractures were modeled explicitly using the EDFM model. To solve the flow equation consistently when the upscaled permeability was a full tensor, a monotone nonlinear two-point flux approximation (NTPFA) [17] method was used so that the comparative results were not distorted by numerical discretization errors.

Description of the Modeling Workflow
The modeling workflow is illustrated in Figure 1. Given the computational grid and DFN, as shown in Figure 1a, a partition length l p is defined to classify all fractures into 2 sets. Fractures whose lengths are smaller than l p are classified as small-scale fractures, whereas all others are classified as large-scale fractures. An illustration is provided in Figure 1b,c. Numerical flow-based upscaling is then used to homogenize the small-scale fractures with the rock matrix, resulting in equivalent permeability tensors for the matrix grid blocks shown in Figure 1d, where the equivalent permeability tensors are represented by ellipses whose semi-axes are scaled by the square root of the principal values of the tensors. The large-scale fractures are, however, modeled explicitly, as shown in Figure 1e. Depending on the local boundary conditions, used for numerical upscaling (homogenization), diagonal or full permeability tensors are obtained. The linear flux approximation scheme, i.e., TPFA, is used to discretize the flow equation if the upscaled tensor is found to be diagonal, and the NTPFA method is used if the upscaled permeability is a full tensor. Note that this leads to an adaptive framework, in which NTPFA is applied only when it is necessary. For both cases, as shown in Figure 1e, the EDFM method is used to model the large-scale fractures along with the heterogeneous matrix rock.

Discrete Fracture Network (DFN) Generation
The DFNs are generated stochastically by drawing samples randomly from specified distributions. Each fracture segment can be characterized by three parameters, namely, its location (the centroid of the line segment c ), orientation ( ), and length ( ). We assume that the distributions of the three parameters are independent. The fracture location c is assumed to follow the uniform distribution. The von Mises distribution is used for fracture orientation . The von Mises distribution is an analogue of the normal distribution for circular data and its probability density distribution function is given by [18] as follows: where ∈ [0, 2 ). Moreover, is the modified Bessel function of the first kind of order zero, and is the distribution mean. The variance of the distribution is given by . The greater the value of , the smaller the variance. The fracture length is often characterized by a positive skewed distribution that has a long tail on the right such as the negative exponential distribution, log-normal distribution, power-law distribution, etc. In this work, the power-law distribution is used to describe the fracture length statistics although other distributions can be used as well. The probability density function for power-law distribution is given by [19] as follows:

Discrete Fracture Network (DFN) Generation
The DFNs are generated stochastically by drawing samples randomly from specified distributions. Each fracture segment can be characterized by three parameters, namely, its location (the centroid of the line segment x c ), orientation (θ), and length (l). We assume that the distributions of the three parameters are independent. The fracture location x c is assumed to follow the uniform distribution. The von Mises distribution is used for fracture orientation θ. The von Mises distribution is an analogue of the normal distribution for circular data and its probability density distribution function is given by [18] as follows: where θ ∈ [0, 2π). Moreover, I 0 is the modified Bessel function of the first kind of order zero, and µ is the distribution mean. The variance of the distribution is given by σ 2 = 1 − I 0 (κ) . The greater the value of κ, the smaller the variance. The fracture length l is often characterized by a positive skewed distribution that has a long tail on the right such as the negative exponential distribution, log-normal distribution, power-law distribution, etc. In this work, the power-law distribution is used to describe the fracture length statistics although other distributions can be used as well. The probability density function for power-law distribution is given by [19] as follows: where l ∈ [l min , l max ]. l min and l max are the minimum and maximum lengths of the fractures, respectively. In addition, α is the exponent that is greater than 1 and C is a normalization constant. Integrating Equation (2) from l min to l max gives the following: Solving Equation (3) for the normalization constant results in C = −α+1 l −α+1 max −l −α+1 min . Note that the stochastically generated DFN is cropped onto a given domain by cutting fractures that intersect the boundary of the domain. Therefore, the actual minimum and maximum fracture lengths are different from the two parameters l min and l max .

Numerical Upscaling
There is a tremendous amount of literature devoted to the subject of upscaling. For the following, we used the techniques presented in [20] and the references therein. Numerically upscaling small-scale fractures with the rock matrix entails solving the incompressible single-phase flow equation on the domain of question. Different upscaling results can be obtained by using different boundary conditions. We considered two types of boundary conditions, namely, constant pressure drop (CP) and linear pressure (LP) boundary conditions. The most common boundary condition is the constant pressure drop on a pair of opposing faces and no-flow boundary conditions on the remaining faces as shown in the left of Figure 2. Suppose that the physical dimension of the grid block to be upscaled is L x , L y in the x, y direction, respectively, for a 2D problem. To compute an equivalent permeability in the x direction, an incompressible single-phase flow problem is solved with a constant pressure drop in the x direction while the remaining boundary faces are closed to flow. Upon computing the total flux Q x in the x direction, the equivalent permeability k * x in the x direction is related to Q x by Darcy's law as follows: where µ is fluid viscosity, A x is the cross-section area for flow in the x direction and is numerically equal to L y in 2D, and ∆p is the pressure drop. Solving for k * x gives the following: In actual implementation, the fluid viscosity µ can be taken as unity and p in = 1, p out = 0 for simplicity. Then k * x is equal to Q x L x /L y . The equivalent permeability k * y in the y direction can be computed similarly by imposing a constant pressure drop in the y direction. Applying the boundary condition of constant pressure drop results in an equivalent permeability that is diagonal. Thus, the classical linear two-point flux approximation (TPFA) method can be used for solving the flow equation for the upscaled model. For brevity, this type of upscaled model is denoted by "Ups-CP". Considering the fact that the discrete fracture network usually has a predominant orientation because of the state of the stress field when the fractures were formed, it may be more appropriate to represent the homogenized media by full permeability tensors with non-zero off-diagonal elements. Therefore, we also applied the linear boundary condition for upscaling, as shown in the right of Figure 2. If the matrix domain and the small-scale fractures contained within it behave as a homogeneous anisotropic media with an equivalent permeability tensor * , then Darcy's equation can be written in a component-wise form as follows: where , are components of the Darcy velocity in the , directions, respectively. * , * , * , * are the components of K*. The total flux in the , direction for single-phase flow is then given by the following: To compute the components of K*, a linear pressure boundary condition ( , ) = 1 − is applied on all the boundary faces of the domain and the incompressible single-phase flow equation is solved. Integrating the Darcy velocity along the respective boundaries gives the total flux , .
Assuming that the pressure distribution inside the domain is also linear and is consistent with the boundary conditions, we can compute the components of the pressure gradient as = − and = 0. Substituting the total flux and pressure gradient into Equation (7) and solving for * , * leads to the following: Considering the fact that the discrete fracture network usually has a predominant orientation because of the state of the stress field when the fractures were formed, it may be more appropriate to represent the homogenized media by full permeability tensors with non-zero off-diagonal elements. Therefore, we also applied the linear boundary condition for upscaling, as shown in the right of Figure 2. If the matrix domain and the small-scale fractures contained within it behave as a homogeneous anisotropic media with an equivalent permeability tensor K * , then Darcy's equation can be written in a component-wise form as follows: where v x , v y are components of the Darcy velocity in the x, y directions, respectively. k * xx , k * xy , k * yx , k * yy are the components of K*. The total flux in the x, y direction for single-phase flow is then given by the following: To compute the components of K*, a linear pressure boundary condition p(x, y) = 1 − x L x is applied on all the boundary faces of the domain and the incompressible single-phase flow equation is solved. Integrating the Darcy velocity along the respective boundaries gives the total flux Q x , Q y . Assuming that the pressure distribution inside the domain is also linear and is consistent with the boundary conditions, we can compute the components of the pressure gradient as ∂p ∂x = − 1 L x and ∂p ∂y = 0. Substituting the total flux and pressure gradient into Equation (7) and solving for k * xx , k * yx leads to the following: where, similar as before, the fluid viscosity is unity. Repeating the same procedure and replacing the boundary condition by p(x, y) = 1 − y L y , we can compute the values of k * xy and k * yy . One issue with this upscaling approach is that the total fluxes Q x , Q y are not uniquely found across the local domain. For example, the total flux computed across the left boundary Q x1 can be different than the one computed across the right boundary Q x2 . For homogeneous isotropic local domains, the two values are close to each other. However, they can differ significantly for heterogeneous anisotropic domains, which is especially true of fractured porous media. In this work, we computed the total flow rates in the x and y directions, i.e., Q x and Q y , as the average of their corresponding end boundaries.
The upscaled full permeability tensors are used to solve the flow equations globally using the nonlinear finite volume method (NTPFA). Derivation of the nonlinear finite volume methods, however, requires that the permeability tensor be symmetric and positive-definite (SPD). The upscaled permeability tensor using the linear pressure boundary condition cannot be guaranteed to be either symmetric or positive-definite. For example, it is not difficult to see that the cross-terms k * xy , k * yx are not equal in general. To enforce symmetry, the final upscaled permeability tensor can be taken as , where the superscript T denotes the transpose. The property of being positive-definite is more difficult to achieve. To circumvent this problem, we checked the positive-definiteness of the upscaled permeability tensor for each upscaled grid block. If it was not positive-definite, we replaced the upscaled permeability with the diagonal tensor by applying the constant pressure drop boundary conditions. Therefore, this second type of upscaled model is actually a hybrid model that utilizes both types of boundary conditions depending on the upscaling results. We refer to this model as "Ups-LP" hereafter.
To solve the incompressible single-phase flow in a fractured porous medium, either the DFM or EDFM models can be used. On the one hand, DFM models require the mesh to conform to the fracture geometries, and unstructured triangular meshes are needed in general. As a result, the mesh may not be K-orthogonal even for isotropic matrix permeability and the linear TPFA method will be inconsistent. Therefore, we used the nonlinear TPFA method for DFM models [21,22], and we denote it by "DFM-NTPFA. EDFM models, on the other hand, work on a Cartesian mesh and the linear TPFA is adequate if the matrix permeability is isotropic. We denote this method by "EDFM-TPFA". The left plot of Figure 3 shows an example of a stochastically generated DFN that has a predominant orientation from the northeast to the southwest on the unit square and the corresponding unstructured triangular mesh is shown in the right of the figure. As can be seen, very fine meshes are needed around fractures to conform to their geometry, especially near fracture intersections. The total number of cells in the mesh is 11,448. Figure 4 shows the pressure distribution of DFM-NTPFA and EDFM-TPFA for the linear pressure boundary condition p(x, y) = 1 − x. The EDFM-TPFA model is solved on a uniform 70 × 70 Cartesian mesh. There is a very close agreement in the pressure solutions between the two methods. Table 1 lists the total flux in the x, y directions computed from the flow solutions of the two methods, and we can see that the results are in excellent agreement. We conducted a series of numerical tests with different realizations of DFNs and found that the results of EDFM-TPFA are quite close to those of DFM-NTPFA for engineering purposes. Therefore, in the following, only the EDFM-TPFA method is used for upscaling as it is computationally more efficient.

Results and Discussion
A series of numerical tests were conducted to assess the performance of the modeling workflow. Whenever possible, the reference solution was obtained by representing all the fractures in the given DFN explicitly using the EDFM method.

Case 1: Effect of Partition Length
In this first test, we established the base case and investigated the effect of partition length on the upscaled model. The reservoir domain was a square with dimensions 500 m × 500 m. The parameters used in this example are listed in Table 2. The left plot of Figure 5 shows one realization of the DFN generated using the given parameters and the right plot of the figure shows the length distribution of the fractures. The actual minimum and maximum fracture lengths were 2.95 m and 360.35 m, respectively, spanning two orders of magnitude. The DFN had a predominant orientation of 45° counterclockwise rotation with respect to the -axis. The total number of fractures in the

Results and Discussion
A series of numerical tests were conducted to assess the performance of the modeling workflow. Whenever possible, the reference solution was obtained by representing all the fractures in the given DFN explicitly using the EDFM method.

Case 1: Effect of Partition Length
In this first test, we established the base case and investigated the effect of partition length on the upscaled model. The reservoir domain was a square with dimensions 500 m × 500 m. The parameters used in this example are listed in Table 2. The left plot of Figure 5 shows one realization of the DFN generated using the given parameters and the right plot of the figure shows the length distribution of the fractures. The actual minimum and maximum fracture lengths were 2.95 m and 360.35 m, respectively, spanning two orders of magnitude. The DFN had a predominant orientation of 45° counterclockwise rotation with respect to the -axis. The total number of fractures in the

Results and Discussion
A series of numerical tests were conducted to assess the performance of the modeling workflow. Whenever possible, the reference solution was obtained by representing all the fractures in the given DFN explicitly using the EDFM method.

Case 1: Effect of Partition Length
In this first test, we established the base case and investigated the effect of partition length on the upscaled model. The reservoir domain was a square with dimensions 500 m × 500 m. The parameters used in this example are listed in Table 2. The left plot of Figure 5 shows one realization of the DFN generated using the given parameters and the right plot of the figure shows the length distribution of the fractures. The actual minimum and maximum fracture lengths were 2.95 m and 360.35 m, respectively, spanning two orders of magnitude. The DFN had a predominant orientation of 45 • counterclockwise rotation with respect to the x-axis. The total number of fractures in the DFN was 2995. All the fractures had the same aperture a, and fracture permeability is given by the following: Dirichlet boundary conditions of 10 MPa were imposed on the bottom left and 5 MPa on the upper right. The reservoir was meshed by a uniform 50 × 50 Cartesian grid. The median of fracture lengths was 7.68 m, meaning that more than half of the fractures had lengths smaller than the characteristic length of the grid. The reference solution was obtained by modeling all the fractures explicitly using the EDFM method. The left plot of Figure 6 shows the reference pressure solution, and the tracer concentration solution after two pore-volume injection (2 PVI) is shown in the right plot of the figure.   The rock matrix was assumed to be homogeneous and isotropic. To drive the global flow, Dirichlet boundary conditions of 10 MPa were imposed on the bottom left and 5 MPa on the upper right. The reservoir was meshed by a uniform 50 × 50 Cartesian grid. The median of fracture lengths was 7.68 m, meaning that more than half of the fractures had lengths smaller than the characteristic length of the grid. The reference solution was obtained by modeling all the fractures explicitly using the EDFM method. The left plot of Figure 6 shows the reference pressure solution, and the tracer concentration solution after two pore-volume injection (2 PVI) is shown in the right plot of the figure.
Energies 2020, 13, x FOR PEER REVIEW 10 of 21 To upscale the small-scale fractures, different partition lengths p were chosen. A smaller value of p means that more fractures are modeled explicitly. Since the upscaling of fractured porous media is tightly related to the size of the domain to be upscaled, we chose the value of p based on the characteristic length of the grid g , as indicated in Table 3. When the value of p is at its largest, only about 1% of fractures are modeled explicitly. Figure 7 shows the relative error of the pressure solutions of the two models, Ups-CP and Ups-LP, for two different partition lengths. It can  To upscale the small-scale fractures, different partition lengths l p were chosen. A smaller value of l p means that more fractures are modeled explicitly. Since the upscaling of fractured porous media is tightly related to the size of the domain to be upscaled, we chose the value of l p based on the characteristic length of the grid l g , as indicated in Table 3. When the value of l p is at its largest, only about 1% of fractures are modeled explicitly. Figure 7 shows the relative error of the pressure solutions of the two models, Ups-CP and Ups-LP, for two different partition lengths. It can be seen that when the partition length is small (l p = l g in the first row), the two models have similar performance in terms of pressure solutions with the maximum relative error less than 1.5%. For a large partition length (l p = 8 × l g in the second row), however, the Ups-LP model performs much better than the Ups-CP model. This can be explained by the fact that when the partition length is smaller than or equal to the characteristic length of the grid, most of the upscaled fractures are contained inside a grid block and their effect on flow is local. Since the rock matrix is homogeneous and isotropic, an upscaled diagonal tensor is adequate to represent the combined flow behavior of the rock matrix and the fractures contained inside it. As the partition length increases, more fractures need to be upscaled with many fractures intersecting domain boundaries. The predominant orientation of the DFN means that the flow has a more pronounced anisotropic behavior at the grid-block scale and a full tensor is more appropriate to represent the fracture media at the scale of interest. Table 3. Partition length for upscaling and the corresponding percentage of fractures modeled explicitly.

Partition Length l p % of Fractures Modeled Explicitly
1 × l g 28.38 2 × l g 9.68 4 × l g 3.37 6 × l g 1.97 8 × l g 1.04 Figure 8 shows the absolute error of the tracer concentration after 2 PVI, and Figure 9 shows the tracer concentration at the outlet as a function of PVI for the two models with the two partition lengths. We can see that, for a small partition length (l p = l g ), the errors of concentration distributions are generally small for both models with the Ups-CP model performing slightly better than the Ups-LP model compared to the reference solution. As the partition length l p increases to 8 × l g , the absolute error of concentration solution also increases for both models as more fractures are homogenized. However, the concentration solution using the Ups-LP model lies closer to the reference solution than the Ups-CP model, which is consistent with observations made about the pressure solutions. Figure 10 further shows the tracer breakthrough curves at the outlet for various partition lengths for the two models. For both models, there is a general tendency for the solutions to converge to the full model solution, i.e., the reference solution, with the Ups-CP model being more sensitive to the partition length. This test case suggests that, when the partition length is small (for example, comparable to the characteristic length of the grid), the Ups-CP model should be used as it is sufficiently accurate and computationally more efficient. For large partition values, the Ups-LP model should be used when the DFN has a predominant orientation as it better captures the directional dependencies of the fracture system. p Explicitly 1 × g 28.38 2 × g 9.68 4 × g 3.37 6 × g 1.97 8 × g 1.04    Ups-LP (second column) models when the partition is l p = l g (first row) and l p = 8 × l g (second row).

Figure 8.
Absolute error of the tracer concentration solutions after 2 PVI for the Ups-CP (first column) and Ups-LP (second column) models when the partition is p = g (first row) and p = 8 × g (second row).

Figure 9.
Tracer concentration at the outlet as a function of pore-volume injected for full model, Ups-CP, and Ups-LP for p = g (left) and p = 8 × g (right).

Case 2: Effect of Fracture Length Distribution
The previous case features a DFN in which more than half of the fractures have lengths smaller than the characteristic length of the grid blocks. For this example, we investigated the effect of fracture length distribution by considering a DFN whose minimum fracture length was larger than the characteristic length of the grid blocks. The parameters used for generating the DFN were the same as those listed in Table 2 with the exception of min . The minimum length here was set to 20 m, which is about twice the size of the characteristic length of the grid. The left plot in Figure 11 shows the generated DFN used in the following simulations, and its histogram of fracture length distribution is shown in the right of the figure. Other parameters for the model setup were the same as the previous case. Figure 12 shows the pressure solution and tracer concentration distribution after 2 PVI for the full model where all the fractures are modeled explicitly using the EDFM method. Unsurprisingly, the system displays a strong preferential orientation for flow and transport because of the presence of the fracture network. For illustration purposes, we first upscaled all the fractures on a 8 × 8 Cartesian grid using the two different upscaling methods, and the results are shown in Figure 13. As expected, Ups-CP produces diagonal tensors, whereas the upscaled tensors from Ups-LP have principal directions aligned with the predominant direction of the underlying DFN. Next, a partition length of p = 5 × g was chosen so that about 10% of the largest fractures were modeled explicitly, whereas the remaining fractures were upscaled with the rock matrix. Figures 14 and 15 display the relative error of the pressure solutions and the absolute error of tracer concentration distribution after 2 PVI for the Ups-CP and Ups-LP models, Figure 10. Tracer concentration at the outlet as a function of pore-volume injected using different partition lengths l p for Ups-CP (left) and Ups-LP (right).

Case 2: Effect of Fracture Length Distribution
The previous case features a DFN in which more than half of the fractures have lengths smaller than the characteristic length of the grid blocks. For this example, we investigated the effect of fracture length distribution by considering a DFN whose minimum fracture length was larger than the characteristic length of the grid blocks. The parameters used for generating the DFN were the same as those listed in Table 2 with the exception of l min . The minimum length here was set to 20 m, which is about twice the size of the characteristic length of the grid. The left plot in Figure 11 shows the generated DFN used in the following simulations, and its histogram of fracture length distribution is shown in the right of the figure. Other parameters for the model setup were the same as the previous case. Figure 12 shows the pressure solution and tracer concentration distribution after 2 PVI for the full model where all the fractures are modeled explicitly using the EDFM method. Unsurprisingly, the system displays a strong preferential orientation for flow and transport because of the presence of the fracture network. For illustration purposes, we first upscaled all the fractures on a 8 × 8 Cartesian grid using the two different upscaling methods, and the results are shown in Figure 13. As expected, Ups-CP produces diagonal tensors, whereas the upscaled tensors from Ups-LP have principal directions aligned with the predominant direction of the underlying DFN. Next, a partition length of l p = 5 × l g was chosen so that about 10% of the largest fractures were modeled explicitly, whereas the remaining fractures were upscaled with the rock matrix. Figures 14 and 15 display the relative error of the pressure solutions and the absolute error of tracer concentration distribution after 2 PVI for the Ups-CP and Ups-LP models, respectively. Figure 14 shows that the relative error of the pressure solutions using the Ups-LP model is significantly smaller than that using the Ups-CP model, suggesting that the Ups-LP model is much more effective in capturing the directional dependencies of the pressure solution caused by the underlying fracture network. The absolute error of the concentration solution shown in Figure 15 show that the Ups-LP model performs modestly better than the Ups-CP model, but the breakthrough curves of tracers at the outlet depicted in Figure 16 shows that the results of Ups-LP are much closer to the reference solution compared to those of the Ups-CP model. The results of this case are in line with the observations made from the first case in that the Ups-LP model is more appropriate for upscaling fractured media when the underlying DFN has a predominant orientation and the upscaled fractures are larger than the characteristic length of the grid. Since the minimum length of the discrete fractures is larger than the characteristic length of the underlying grid, the fractures to be upscaled generally cut through the hosting grid block, which, combined with the fact that the fracture network has a predominant orientation, means that the Ups-LP model has superior performance compared to the Ups-CP model.

Case 3: Effect of Fracture Orientation
For this case, we investigated the effect of fracture orientation on the performance of the two upscaling methods. The parameters used for generating the DFN are the same as those listed in Table 2 except that the parameter takes the value of 1. As mentioned in Section 2.2, the parameter controls the variance of the fracture orientations. A small value of 1 gives an orientation distribution that has large variance. Figure 17 shows the corresponding DFN used in the following simulations together with the histogram of the length distribution. As can be seen from the figure, the DFN does not have a particular predominant orientation. Similarly, the reference solution is obtained by modeling all the fractures explicitly using the EDFM method, and the results are given in Figure 18, which shows that the pressure diffuses more uniformly from the injection inlet to the production outlet and the distribution of tracer concentration is also more random compared to the previous two cases. Figures 19 and 20 show the relative errors of pressure solutions and absolute errors of tracer concentration solutions after 2 PVI, respectively, for the two upscaling models, Ups-CP and Ups-LP. For this case, the Ups-CP model shows better performance in terms of pressure solutions over the Ups-LP model as can be seen in Figure 19, but conclusions about the transport solutions are less definitive from Figure 20. Furthermore, there is no obvious winner that performs significantly better than the other, although the breakthrough curves depicted in Figure 21 suggest that the flux solution of Ups-LP is closer to the reference solution than that of Ups-CP. Different realizations of the DFN were generated using the same parameters and additional simulations were carried out on these DFNs. Our findings indicated that the opposite scenario can also occur; that is, for some realizations of the DFN, Ups-LP gives more accurate pressure solutions but less accurate flux solutions than Ups-CP compared to the reference solution. Therefore, no definite conclusions can be reached about which one is better when the DFN does not have predominant orientations.

Case 3: Effect of Fracture Orientation
For this case, we investigated the effect of fracture orientation on the performance of the two upscaling methods. The parameters used for generating the DFN are the same as those listed in Table 2 except that the parameter κ takes the value of 1. As mentioned in Section 2.2, the parameter κ controls the variance of the fracture orientations. A small value of 1 gives an orientation distribution that has large variance. Figure 17 shows the corresponding DFN used in the following simulations together with the histogram of the length distribution. As can be seen from the figure, the DFN does not have a particular predominant orientation. Similarly, the reference solution is obtained by modeling all the fractures explicitly using the EDFM method, and the results are given in Figure 18, which shows that the pressure diffuses more uniformly from the injection inlet to the production outlet and the distribution of tracer concentration is also more random compared to the previous two cases. Figures 19 and 20 show the relative errors of pressure solutions and absolute errors of tracer concentration solutions after 2 PVI, respectively, for the two upscaling models, Ups-CP and Ups-LP. For this case, the Ups-CP model shows better performance in terms of pressure solutions over the Ups-LP model as can be seen in Figure 19, but conclusions about the transport solutions are less definitive from Figure 20. Furthermore, there is no obvious winner that performs significantly better than the other, although the breakthrough curves depicted in Figure 21 suggest that the flux solution of Ups-LP is closer to the reference solution than that of Ups-CP. Different realizations of the DFN were generated using the same parameters and additional simulations were carried out on these DFNs. Our findings indicated that the opposite scenario can also occur; that is, for some realizations of the DFN, Ups-LP gives more accurate pressure solutions but less accurate flux solutions than Ups-CP compared to the reference solution. Therefore, no definite conclusions can be reached about which one is better when the DFN does not have predominant orientations.

Summary and Conclusions
A general workflow for modeling flow and transport in fractured porous media was introduced in this work. Based on a given partition length p, fractures were classified as small-scale fractures (fractures with lengths smaller than p) and large-scale fractures (fractures with lengths greater than p). The small-scale fractures were upscaled with the hosting matrix, whereas the large-scale fractures were modeled explicitly using the EDFM model. A hybrid upscaling approach was proposed that switches between applying the linearly varying pressure boundary condition and the constant pressure drop boundary condition. The linearly varying pressure boundary condition was preferentially applied to capture the anisotropic effect of the fracture porous media.

Summary and Conclusions
A general workflow for modeling flow and transport in fractured porous media was introduced in this work. Based on a given partition length l p , fractures were classified as small-scale fractures (fractures with lengths smaller than l p ) and large-scale fractures (fractures with lengths greater than l p ).
The small-scale fractures were upscaled with the hosting matrix, whereas the large-scale fractures were modeled explicitly using the EDFM model. A hybrid upscaling approach was proposed that switches between applying the linearly varying pressure boundary condition and the constant pressure drop boundary condition. The linearly varying pressure boundary condition was preferentially applied to capture the anisotropic effect of the fracture porous media. If the resulting upscaled permeability tensor was not positive-definite, the constant pressure drop boundary condition was used instead to produce an upscaled tensor that was diagonal. Once the small-scale fractures were upscaled, the original fractured reservoir was replaced by an equivalent medium with large-scale fractures embedded in it. Our numerical results suggest that when the partition length is comparable to the characteristic length of the grid blocks or the DFN does not have a predominant orientation, the Ups-CP model should be the method of choice as it is computationally more efficient and has similar accuracy compared to the Ups-LP model. The Ups-LP model is best suited for cases in which the small-scale fractures to be upscaled cut through the grid block and have a predominant orientation.
The current workflow implementation homogenizes the small-scale fractures and the porous matrix into a single medium. An alternative is to upscale only the small-scale fractures into a distinct continuum of its own, leading to dual-continuum models while the large-scale fractures are still modeled explicitly. A comparison between the single-continuum upscaling and the dual-continuum upscaling would be of great interest for obtaining a complete picture of modeling small-scale fractures. Moreover, the current work focuses on single-phase flow. A future extension could include multi-phase flow along with its intertwined complex flow dynamics.