Discrete Fracture-matrix Model Research Articles

Quantifying the impact of upscaled parameters on radionuclide transport in three-dimensional fracture-matrix systems

Assessing the long-term safety of geological repositories for high-level radioactive waste is critically dependent on understanding radionuclide transport in multi-scale fractured rocks. This study explores the influence of upscaled parameters on radionuclide movement within a three-dimensional fracture-matrix system using a discrete fracture-matrix (DFM) model. The developed numerical simulation workflow includes creating a random discrete fracture network, meshing of the fractures and matrix, assigning upscaled parameters, and conducting finite element simulations. We simulated the spatiotemporal evolution of radionuclide concentrations in the fractures and matrix over a century, revealing significant spatial heterogeneity driven by a heterogeneous seepage field. Employing geostatistics-based upscaling methods, we predicted the effective ranges of crucial solute transport parameters at the field scale. The matrix diffusion coefficient, matrix distribution coefficient, and longitudinal dispersivity were upscaled by factors of 2.0–3.0, 2.5–4.0, and 10–104, respectively, based on laboratory-scale measurements. Incorporating these upscaled parameters into the DFM model, we analyzed their impact on radionuclide transport. Our findings demonstrate that an upscaled matrix diffusion coefficient and matrix distribution coefficient result in a delayed transport of radionuclides in fractures by enhancing mass transfer between the fractures and rock matrix, while an upscaled longitudinal dispersivity accelerates transport by advancing the positions of concentration peaks. Sensitivity analysis revealed that the matrix distribution coefficient is the most impactful, followed by dispersivity and matrix diffusion coefficient. These insights are important for minimizing parameter uncertainties and enhancing the accuracy of predictions concerning radionuclide transport in multi-scale fractured rocks.

Simulating multiphase flow in fractured media with graph neural networks

Numerical modeling of flow dynamics with multiple fluid phases in subsurface fractured porous media is of great significance to numerous geoscience applications. Discrete fracture-matrix (DFM) approaches become popular for simulating fractured reservoirs in the last decade. Data-driven surrogate models can provide computationally efficient alternatives to high-fidelity numerical simulators. Although convolutional neural networks (CNNs) are effective at approximating the space-time solutions of multiphase flowing processes, it remains difficult for CNNs to operate upon DFMs with unstructured meshes. To tackle this challenge, we leverage graph neural networks (GNNs) for surrogate modeling of an embedded DFM model. The results of two-dimensional cases with complex fracture systems show that the learned surrogates precisely capture the effect of the variations in fracture connectivity and forecast dynamic pressure and saturation solutions with high accuracy. Furthermore, our GNN-based models exhibit promising generalizability to fracture networks with different geometries and numbers of fractures that are not encountered from the training dataset.

Numerical simulation of CO2-enhanced oil recovery in fractured shale reservoirs using discontinuous and continuous Galerkin finite element methods

Introduction: This study explores the potential of enhancing shale oil recovery and reducing CO2 emissions through CO2 injection in fractured shale reservoirs. The importance of this approach lies in its dual benefit: improving oil extraction efficiency and addressing environmental concerns associated with CO2 emissions.Method: We employed a discrete fracture-matrix model to simulate CO2 flooding in fractured shale reservoirs, utilizing both discontinuous Galerkin (DG) and continuous Galerkin (CG) finite element methods. The DG-CG FEM’s accuracy was validated against the McWhorter problem, ensuring the reliability of the simulation results. Our model also considered various factors, including reservoir heterogeneity, fracture permeability, CO2 injection volume, and gas injection patterns, to analyze their impact on shale oil recovery.Result: Our simulations revealed that fractured reservoirs significantly enhance shale oil production efficiency compared to homogeneous reservoirs, with an approximate 48.9% increase in production. A notable increase in shale oil production, by 15.8%, was observed when fracture permeability was increased by two orders of magnitude. Additionally, a fourfold increase in CO2 injection rate resulted in a 31.5% rise in shale oil production. Implementing a step-by-step reduction in injection volume while maintaining the total CO2 injection constant proved to be more effective than constant-rate injections.Discussion: The study demonstrates the effectiveness of CO2 flooding in fractured shale reservoirs for enhancing shale oil recovery.

An efficient 3D cell-based discrete fracture-matrix flow model for digitally captured fracture networks.

Complex hydraulic fracture networks are critical for enhancing permeability in unconventional reservoirs and mining industries. However, accurately simulating the fluid flow in realistic fracture networks (compared to the statistical fracture networks) is still challenging due to the fracture complexity and computational burden. This work proposes a simple yet efficient numerical framework for the flow simulation in fractured porous media obtained by 3D high-resolution images, aiming at both computational accuracy and efficiency. The fractured rock with complex fracture geometries is numerically constructed with a cell-based discrete fracture-matrix model (DFM) having implicit fracture apertures. The flow in the complex fractured porous media (including matrix flow, fracture flow, as well as exchange flow) is simulated with a pipe-based cell-centered finite volume method. The performance of this model is validated against analytical/numerical solutions. Then a lab-scale true triaxial hydraulically fractured shale sample is reconstructed, and the fluid flow in this realistic fracture network is simulated. Results suggest that the proposed method achieves a good balance between computational efficiency and accuracy. The complex fracture networks control the fluid flow process, and the opened natural fractures behave as primary fluid pathways. Heterogeneous and anisotropic features of fluid flow are well captured with the present model.

Impact of artificial topological changes on flow and transport through fractured media due to mesh resolution

We performed a set of numerical simulations to characterize the interplay of fracture network topology, upscaling, and mesh refinement on flow and transport properties in fractured porous media. We generated a set of generic three-dimensional discrete fracture networks at various densities, where the radii of the fractures were sampled from a truncated power-law distribution, and whose parameters were loosely based on field site characterizations. We also considered five network densities, which were defined using a dimensionless version of density based on percolation theory. Once the networks were generated, we upscaled them into a single continuum model using the upscaled discrete fracture matrix model presented by Sweeney et al. (2019). We considered steady, isothermal pressure-driven flow through each domain and then simulated conservative, decaying, and adsorbing tracers using a pulse injection into the domain. For each simulation, we calculated the effective permeability and solute breakthrough curves as quantities of interest to compare between network realizations. We found that selecting a mesh resolution such that the global topology of the upscaled mesh matches the fracture network is essential. If the upscaled mesh has a connected pathway of fracture (higher permeability) cells but the fracture network does not, then the estimates for effective permeability and solute breakthrough will be incorrect. False connections cannot be eliminated entirely, but they can be managed by choosing appropriate mesh resolution and refinement for a given network. Adopting octree meshing to obtain sufficient levels of refinement leads to fewer computational cells (up to a 90% reduction in overall cell count) when compared to using a uniform resolution grid and can result in a more accurate continuum representation of the true fracture network.

DFNMesh: Finite element meshing for discrete fracture matrix models

This paper presents a novel methodology for robust Finite Element (FE) mesh generation of Discrete Fracture Matrix (DFM) models. The method can handle the complete multi-dimensional domain from the tridimensional porous rock matrix, through fracture surfaces, down to open curves for fracture–fracture intersections. The accompanying open-source code is written in modern C++ with a JSON interface and largely relies on two state-of-the-art FE libraries: NeoPZ and Gmsh. Starting from a user-defined coarse mesh, fractures are sequentially read as convex polygons and introduced without distortions to the coarse geometry. The main steps involve: intersect edges by checking for nodes on opposite sides of the fracture plane, extend intersections from edges to faces, coalesce intersections to closest existing nodes (given a tolerance), refine interface elements to conform to the fracture, identify subsets of fracture surface, mesh the surface, and locate boundaries and intersections where they arise. Finally, the space around fractures is filled with the fine-scale unstructured mesh, which is kept conformal. The robustness of the implementation is derived from the consistent background of well-defined and simple premises like convexity and side-specific element neighborhood. Results show that the proposed technique can construct adequate 3D DFM grids, while still giving users freedom to adjust between geometrical fidelity and mesh quality through more aggressive feature rejection.

Computed-tomography-based discrete fracture-matrix modeling: An integrated framework for deriving fracture networks

Investigating the topology of natural fractures and associated flow properties allow for better understanding of fluid flow patterns and mechanical response in fractured reservoirs, key points in the energy supply through geothermal systems, or in the reduction of net CO2 emissions through carbon capture and storage (CCS) practices. However, limited information on the fracture networks is available from seismic data, which poses as a challenge in the accurate modeling of fractured reservoirs. This work proposes an integrated framework for obtaining fracture networks with real aperture distributions, based on X-ray computed tomography (CT) images of core plugs from naturally fractured rock formations. We adapt and employ segmentation, skeletonization and geometry generation algorithms that are available in open-source libraries and licensed software. In addition, we investigate the contributions of fractures and porous matrix to the equivalent permeability tensor through fluid flow simulations considering three different scenarios for the matrix permeability. This work is novel in the coordination of the steps involved in the generation of discrete fracture-matrix (DFM) models from CT-scans of core plugs.

Efficient dissipation-based nonlinear solver for multiphase flow in discrete fractured media

Solving nonlinear equation systems resulting from the Fully Implicit Method (FIM) remains a challenge for numerical simulation of multiphase flow in subsurface fractured media. The Courant numbers can vary orders of magnitude across discrete fracture-matrix (DFM) models because of the high contrasts in the permeability and length-scale between matrix and fracture. The standard Newton solver is usually unable to converge for big timestep sizes.Limited research has been conducted on nonlinear solver techniques for multiphase flow and transport in fractured media. We make an extension of a new dissipation-based continuation (DBC) method to DFM models. Our goal is to prevent timestep cuts and maintain efficient time-stepping for FIM. The DBC algorithm builds a homotopy of the discretized conservation equations through the addition of numerical dissipation terms. We introduce a continuation parameter for controlling the dissipation and ensuring that accuracy of the computed solution will not be reduced. Under the nonlinear framework of DBC, general dissipation operators and adaptive strategies are developed to provide the optimal dissipation matrix for the multiphase hyperbolic system.We assess the new nonlinear solver using multiple numerical examples. Results reveal that the damped-Newton solver suffers from serious restrictions on timestep sizes and wasted iterations. In contrast, the DBC solver provides excellent computational performance. The dissipation operators are able to successfully resolve main convergence difficulties. We also investigate the impact of star-delta transformation which removes the small cells at fracture intersections. Moreover, we demonstrate that an aggressive time-stepping does not affect the solution accuracy.

A robust fully mixed finite element model for flow and transport in unsaturated fractured porous media

A fully mixed finite element (MFE) model is developed for nonlinear flow and transport in unsaturated fractured porous media with matrix-fracture and fracture-fracture fluid and mass exchanges. The model is based on the discrete fracture matrix (DFM) approach and assumes cross-flow equilibrium in the fractures. The MFE method is employed for the spatial discretization of both flow and transport on the 2D-matrix elements as well as on the 1D-fracture elements. An upwind scheme is employed to avoid unphysical oscillations in the case of advection dominant transport. The temporal discretization is performed using high-order time integration methods and efficient automatic time-stepping schemes via the MOL.Two test problems dealing with flow and mass transport in saturated and unsaturated fractured porous media are simulated to show the validity of the new model by comparison against (i) a 1D-2D Comsol finite element model and (ii) a 2D-2D Discontinuous Galerkin (DG) model where both fractures and matrix continua are discretized with small 2D mesh elements. The robustness and efficiency of the developed 1D-2D MFE model are then investigated for a challenging problem dealing with infiltration of contaminated water into an initially dry soil involving a fracture network.The new model yields stable results for advection-dominated and advection-dispersion transport configurations. Further, the results of the 1D-2D MFE model are in very good agreement with those of the 2D-2D DG model for both configurations. The simulation of infiltration of contaminated water into a dry fractured soil shows that the 1D-2D MFE model is within 15 times more efficient than the 2D-2D DG model, which confirms the high benefit of using robust and efficient DFM models for the simulation of flow and transport in fractured porous media.

Flow and transport in three-dimensional discrete fracture matrix models using mimetic finite difference on a conforming multi-dimensional mesh

We present a comprehensive workflow to simulate single-phase flow and transport in fractured porous media using the discrete fracture matrix approach. The workflow has three primary parts: (1) a method for conforming mesh generation of and around a three-dimensional fracture network, (2) the discretization of the governing equations using a second-order mimetic finite difference method, and (3) implementation of numerical methods for high-performance computing environments. A method to create a conforming Delaunay tetrahedralization of the volume surrounding the fracture network, where the triangular cells of the fracture mesh are faces in the volume mesh, that addresses pathological cases which commonly arise and degrade mesh quality is also provided. Our open-source subsurface simulator uses a hierarchy of process kernels (one kernel per physical process) that allows for both strong and weak coupling of the fracture and matrix domains. We provide verification tests based on analytic solutions for flow and transport, as well as numerical convergence. We also provide multiple expositions of the method in complex fracture networks. In the first example, we demonstrate that the method is robust by considering two scenarios where the fracture network acts as a barrier to flow, as the primary pathway, or offers the same resistance as the surrounding matrix. In the second test, flow and transport through a three-dimensional stochastically generated network containing 257 fractures is presented.

A finite-volume moving-mesh method for two-phase flow in dynamically fracturing porous media

Multiphase flow in fractured porous media can be described by discrete fracture matrix models that represent the fractures as dimensionally reduced manifolds embedded in the bulk porous medium. Generalizing earlier work on this approach we focus on immiscible two-phase flow in time-dependent fracture geometries, i.e., the fracture itself and the aperture of the fractures might evolve in time. For dynamic fracture geometries of that kind, neglecting capillary forces, we deduce by transversal averaging of a full dimensional description a dimensionally reduced model that governs the geometric evolution and the flow dynamics. The core computational contribution is a mixed-dimensional finite-volume discretization based on a conforming moving-mesh ansatz. This finite-volume moving-mesh (FVMM) algorithm is tracking the fractures' motions as a family of unions of facets of the mesh. Notably, the method permits arbitrary movement of facets of the triangulation while keeping the mass conservation constraint. In a series of numerical examples we investigate the modeling error of the reduced model as it compares to the original full dimensional model. Moreover, we show the performance of the finite-volume moving-mesh algorithm for the complex wave pattern that is induced by the interaction of saturation fronts and evolving fractures.

Matrix acidization in fractured porous media with the continuum fracture model and thermal Darcy-Brinkman-Forchheimer framework

The numerical simulation of matrix acidization in fractured porous media is continuous work in the reservoir simulation community. However, existing works mainly depend on two kinds of models to describe fractured porous media: the discrete fracture-matrix (DFM) model and the continuum fracture (CF) model. In addition, most of the works use Darcy's equation to describe the flows. The DFM model can accurately describe the fracture and the flow in the fracture at the cost of the complexity of the simulation framework. On the other hand, the CF model treats the fracture and the matrix uniformly, which reduces the complexity of the simulation framework, but the accuracy of the results is not as good as that of the DFM model. In summary, the challenge of this area is that the accuracy of the results cannot be guaranteed by a sufficiently simple model. The Darcy-Brinkman-Forchheimer (DBF) framework was originally developed to handle an issue in matrix acidization simulation, i.e., the flows in the high-porosity zones cannot be well described by Darcy's law. Considering the similarity of the fracture and the high-porosity zones, it is anticipated that the proposed DBF framework can also describe the flows in the fracture well., i.e., the accuracy of the results guaranteed by the DFM model can also be guaranteed by the DBF framework. Compared with the DFM model, the DBF framework is much simpler. Moreover, the DBF framework can be naturally integrated with the CF model. Thus, this work leverages the CF model and the DBF framework to simulate matrix acidization in fractured porous media, by which accurate results can be achieved with a simple simulation framework. The simulation results in this work are compared with those of Khoei's work, by which the reasonability of this work is demonstrated. After that, the thermal DBF framework, which is an expansion of the DBF framework, is leveraged to investigate matrix acidization in fractured porous media under different thermal conditions, and many helpful conclusions that benefit realistic operations are drawn. For example, the fracture orientation can have different effects on matrix acidization under different temperature conditions.

Performance Analysis of Multi-Task Deep Learning Models for Flux Regression in Discrete Fracture Networks

In this work, we investigate the sensitivity of a family of multi-task Deep Neural Networks (DNN) trained to predict fluxes through given Discrete Fracture Networks (DFNs), stochastically varying the fracture transmissivities. In particular, detailed performance and reliability analyses of more than two hundred Neural Networks (NN) are performed, training the models on sets of an increasing number of numerical simulations made on several DFNs with two fixed geometries (158 fractures and 385 fractures) and different transmissibility configurations. A quantitative evaluation of the trained NN predictions is proposed, and rules fitting the observed behavior are provided to predict the number of training simulations that are required for a given accuracy with respect to the variability in the stochastic distribution of the fracture transmissivities. A rule for estimating the cardinality of the training dataset for different configurations is proposed. From the analysis performed, an interesting regularity of the NN behaviors is observed, despite the stochasticity that imbues the whole training process. The proposed approach can be relevant for the use of deep learning models as model reduction methods in the framework of uncertainty quantification analysis for fracture networks and can be extended to similar geological problems (for example, to the more complex discrete fracture matrix models). The results of this study have the potential to grant concrete advantages to real underground flow characterization problems, making computational costs less expensive through the use of NNs.

Laboratory and numerical investigation of saline intrusion in fractured coastal aquifers

Laboratory scale experiments and numerical modelling were employed in this study to investigate saltwater intrusion in fractured aquifers. Saline intrusion was initiated in one homogeneous and six fractured experimental aquifers containing individual discontinuities of varying length and orientation. Automated image processing enabled high precision quantification of three intrusion variables, the toe length of the saline wedge, the width of the mixing zone and the aquifer fraction occupied by saltwater. A discrete fracture matrix model was successfully utilized to recreate the experimental data and expand the study's findings through rigorous sensitivity analysis. The presence of fractures significantly impacted all three intrusion variables under consideration. The length of intrusion was negatively correlated to the horizontal fracture's distance from the systems’ seaward boundary. It was demonstrated that for the same fractured aquifer, the presence of a discontinuity can either limit or augment saline intrusion, depending on the applied hydraulic gradient. For gradients steeper than a critical head difference, at which the toe length was the same for both the fractured and the equivalent homogeneous aquifer, intrusion was suppressed further seaward, while for milder ones it intensified. The distance of horizontal fracture from the aquifer's base determined the extend of intrusion in the vertical direction. In general, the longer the discontinuities were, the more significant their impact on groundwater dynamics. In the case of vertical factures, whenever the saline wedges reached their position, the discontinuities contributed significantly in the widening of the mixing zone, while having limited effect on the other two intrusion characteristics. In aquifers with discontinuities adjacent to their side boundaries, a distinct distribution of saltwater concentration was identified, distinguishing them from the rest of the aquifers.

Modeling and Analysis of the Coupling in Discrete Fracture Matrix Models

This paper deals with the derivation and analysis of reduced order elliptic PDE models on fractured domains. We use a Fourier analysis to obtain coupling conditions between subdomains when the frac...

An advanced discrete fracture model for variably saturated flow in fractured porous media

Accurate modeling of variably saturated flow (VSF) in fractured porous media with the discrete fracture-matrix (DFM) model is a computationally challenging problem. The applicability of DFM model to VSF in real field studies at large space and time scales is often limited, not only because it requires detailed fracture characterization, but also as it involves excessive computational efforts. We develop an efficient numerical scheme to solve the Richards equation in discretely fractured porous media. This scheme combines the mixed hybrid finite element method for space discretization with the method of lines for time integration. The fractures are modeled as lower-dimensional interfaces (1D), within the 2D porous matrix. We develop a new mass-lumping (ML) technique for the fractures to eliminate unphysical oscillations and convergence issues in the solution, which significantly improves efficiency, enabling larger field applications. The proposed new scheme is validated against a commercial simulator for problems involving water table recharge at the laboratory scale. The computational efficiency of the developed scheme is examined on a challenging problem for water infiltration in fractured dry soil, and compared with standard numerical techniques. We show that the ML technique is crucial to improve robustness and efficiency, which outperforms the commonly used methods that we tested. The applicability of our method is then demonstrated in a study concerning the effect of climate change on groundwater resources in a karst aquifer/spring system in El Assal, Lebanon. Simulations, including recharge predictions under climate change scenarios, are carried out for about 80 years, up to 2099. This study demonstrates the applicability of our proposed scheme to deal with real field cases involving large time and space scales with high variable recharge. Our results indicate that the water-table level is sensitive to the presence of fractures, where neglecting fractures leads to an overestimation of the available groundwater amount. The proposed numerical approach is generic for DFM model and can be extended to different 2D and 3D finite-element frameworks.

A Discrete fracture–matrix model for pressure transient analysis in multistage fractured horizontal wells with discretely distributed natural fractures

Analysis of transient pressure test has become a reliable and powerful tool for petroleum engineers to study the impacts of complex fracture networks on multistage fractured horizontal wells (MFHWs) in unconventional reservoirs. Although lots of models have been established in the past decades, such as well-known dual-porosity model (DPM), dual-permeability model (DKM) and stripe fracture model (SFM), they have been proven unsuitable for modeling discretely distributed natural fractures (NFs).In this study, to overcome the disadvantages of existing models, a Discrete fracture–matrix method based on numerical well testing model (DFM-WTM) is proposed to study the pressure transient responses of MFHW with discretely distributed NFs. In DFM-WTM, local grid refinement method, well model and transmissibility corrections, natural fracture generation and time-step size are improved based on the framework of Matlab Reservoir Simulation Toolbox (MRST) (Lie, 2019), ensuring that pressure transient responses of MFHW with complex fracture networks can be obtained accurately and efficiently. Then, sensitivity analysis is conducted to study the impacts of different fracture properties. Finally, a case study is performed to further demonstrate the reliability and practicality of DFM-WTM.The results show that MFHW with discretely distributed NFs can be divided into 8 flow regimes, including (1) wellbore storage and skin effect, (2) fracture linear flow, (3) bilinear flow, (4) NF–HF fluid supply, (5) second bilinear flow, (6) formation linear flow, (7) pseudo boundary flow and (8) pseudo radial flow. Besides, a “dip” is found in the pressure derivative curve which is caused by the NF–HF fluid supply. Through the sensitivity analysis, it is found that this “dip” is an important signal to identify the impacts of NFs. It can also help us fathom the properties of HFs and NFs because this “dip” can be mainly affected by the number and conductivity of both HFs and NFs.

Discrete Fracture Matrix Model Applied to the Computation of Water Flow Through the Underground Facility

Mathematical model of groundwater flow in fractured porous media implemented in the GeRa code is presented. It uses the discrete fracture matrix approach which allows for explicit fractures approximation and couples 2D and 3D flow models in fractures and porous matrix respectively. Discretization is done using finite volume method on general conformal polyhedral grids, while 2D surface grids on fractures being the trace of 3D grid in the porous media on the fractures’ surfaces. The model is applied to assess the effect of fracture zones on groundwater flow through a projected radioactive waste deep geological disposal.

Upscaled discrete fracture matrix model (UDFM): an octree-refined continuum representation of fractured porous media

We develop a computational geometry-based upscaling approach to accurately capture the dynamic processes occurring within and between subsurface fracture networks and the surrounding porous matrix named the Upscaled Discrete Fracture Matrix model (UDFM). Fracture attributes (orientations and apertures) are upscaled and combined with matrix attributes (permeability and porosity) into properties of an unstructured spatially variable Delaunay tetrahedral continuum mesh. The resolution of the mesh depends on proximity to the fractures to preserve the geometric and topological integrity of the underlying fracture network as well as increasing the accuracy of gradients in the dynamics between the fractures and matrix. The UDFM Delaunay mesh and its dual Voronoi mesh can be used in existing multiphysics simulators for flow, solute/heat mass transport, and geomechanics, thereby eliminating the need for the additional development of numerical methods to couple processes in fracture/matrix systems. The model is verified by comparing flow, transport, and coupled heat-mass flow simulations on the provided meshes against analytical and numerical benchmarks. We also provide an additional example to demonstrate the applicability of the UDFM to complex heteregenous fractured media. Overall, the UDFM is accurate in all of the cases considered here and presents an attractive alternative to other modeling strategies that require novel numerical methods or multi-dimensional meshes.

Representation of a new physics-based non-Darcy equation for low-velocity flow in tight reservoirs

There has been a growing interest in the study of low-velocity non-Darcy flow in tight reservoirs. Most existing non-Darcy equations, however, are flawed in model assumption, equation derivation or field practicability. For another, current non-Darcy simulations are based on the traditional methods, which are incapable of describing the fracture-rich feature of tight reservoir geology. In light of this, this paper introduces a new physics-based non-Darcy equation for low-velocity flow from four aspects: equation derivation, experimental validation, flow simulation and flow analysis. The rigorous derivations of the new equation are presented with a modified Buckingham-Reiner equation in which the boundary layer effect and the Bingham fluid behavior are described. The non-Darcy coefficient controlled by the boundary layer parameter and threshold pressure gradient (TPG) is incorporated into Darcy's equation. The equation validations are performed by liquid column experiment and core flow experiment, which yields good matching relationship between modeled values and observed values. The non-Darcy flow simulation is based on the triangle-gridded discrete fracture-matrix model (DFM) which can model arbitrary fractures. The upstream weighting scheme for computing the non-Darcy transmissibility is proposed to improve numerical robustness. Then the scheme is implemented by constructing non-Darcy operators, consisting of half-face upstream weighting operator and gradient operator, and further setting up the Automatic Differentiation (AD)-style discrete equations. The studies on non-Darcy flow in matrix coupled with Darcy flow in fracture are performed by modeling two types of fractured reservoirs. Detailed analyses are illustrated with pressure transient, long-term well bottom pressure and pressure map.