Discrete Fracture Model Research Articles

Multiscale model reduction for the time fractional thermoporoelasticity problem in fractured and heterogeneous media

In this paper, we consider the time fractional thermoporoelasticity problem in fractured and heterogeneous media. The mathematical model with a time memory formalism is described by a coupled system of equations for pressure, temperature and displacements. We use an implicit finite difference approximation for temporal discretization. We present a fine grid approximation based on the finite element method and Discrete Fracture Model (DFM) for two-dimensional model problems. Further, we use the Generalized Multiscale Finite Element Method (GMsFEM) for coarse grid approximation. The primary concept behind the proposed method is to streamline the complexity inherent in the thermoporoelasticity problem. Given that our model equation incorporates multiple fractional powers, leading to multiple unknowns with memory effects, we aim to address this intricacy by optimizing the problem’s dimensionality. As a result, the solution is sought on a coarse grid, a strategic choice that not only simplifies the computational cost but also contributes to significant time savings. We present numerical results for the two-dimensional model problems in heterogeneous fractured porous media. We derive relative errors between the reference fine grid solution and the multiscale solution for different numbers of multiscale basis functions. The results confirm that the proposed method is able to achieve good accuracy with a few degrees of freedoms on the coarse grid.

Advancements and Perspectives in Embedded Discrete Fracture Models (EDFM)

The Embedded Discrete Fracture Model (EDFM) has emerged as a prominent piece of technology used for embedding the hydraulic behavior of rock joints in reservoir numerical models. This paper critically reviews its fundamentals, the latest developments, and opportunities for further research. The literature is extensive regarding novel algorithms attempting to reach more accurate and computationally effective estimates. While hydraulic fracture models seem suitable for their purposes, their assumptions might be excessively simplistic and unrealistic when assessing naturally fractured reservoirs. The paper begins by examining fractures as physical characteristics and the key mechanisms to be considered when integrating them into numerical flow simulators. The use of the EDFM technique shows promise for simulating capillary continuity and buoyancy effects in multiphase and multicomponent cases. However, there are significant limitations that hinder its widespread field-scale adoption for reservoir performance evaluation. In this regard, the lack of public-domain realistic benchmarks to validate and compare the potential of each method reinforces the difficulties of performing broader applications of the EDFM techniques in large-scale models.

Numerical simulation on reservoir stimulation assisted depressurization development of hydrate bearing layers based on embedded discrete fractures

As the global demand for energy escalates, the increasing severity of environmental pollution underscores the urgent need for efficient development of hydrates—a promising, clean energy source with vast reserves. It is imperative to carry out reservoir stimulation on low-permeability hydrate reservoirs to enhance fluidity. In order to simulate the geometric morphology of complex fractures, based on the bifurcation fracture morphology in existing laboratory experiments, combined with embedded discrete fracture model and tree fractal theory, the multiphase and multi-component seepage mechanism, thermal conduction, thermal convection, and hydrate phase change in hydrate-bearing layers were comprehensively considered. Based on the TOUGH code, a TOUGH-EDFM module is developed to realize numerical simulation of three-dimensional four-phase four-component hydrate reservoir stimulation assisted depressurization development. We clarify the production capacity change characteristics, main methane dissociation zones, and gas-water migration trends after reservoir stimulation. The results indicate that in the multi-level branching fracture area, the increase in hydrate production in the primary fracture area of the reservoir is greater than that in the secondary fracture branching area. When the fractal scale factor is 0.56, the increase in production in the secondary branch fracture zone is 48% of that in the primary fracture zone. The cumulative gas production after the reservoir stimulation of Class III hydrate reservoirs can surpass that of the basic depressurization scheme by a remarkable 2.21 times, and the overall hydrate dissociation proportion has increased by 60%, indicating the ability of the reservoir to effectively utilize the energy within hydrate reservoirs, thereby significantly improving its depressurization development potential.

A compositional numerical study of vapor–liquid-adsorbed three-phase equilibrium calculation in a hydraulically fractured shale oil reservoir

The development of carbon capture, utilization, and storage technologies has notably advanced CO2-enhanced oil recovery (EOR) in shale oil reservoirs, which are characterized by abundant nanopores. These nanopores induce unique phase behaviors in hydrocarbons, challenging traditional phase equilibrium calculation methods. This paper presents a novel three-phase thermodynamic model (vapor–liquid-adsorbed three-phase equilibrium calculation) that addresses these challenges by considering the nanopore capillary pressure, critical parameter transitions, and material exchange between the adsorbed and bulk phases. Grounded in the multicomponent Langmuir–Freundlich adsorption equation and the Peng Robinson equation of state, this model is integrated into the MATLAB Reservoir Simulation Toolbox using an embedded discrete fracture model framework, enabling detailed study of CO2 and hydrocarbon phase behaviors within shale oil nanopores. The results reveal that there are significant nano-constrained effects on multicomponent fluid phase behavior, particularly in pores smaller than 20 nm, leading to notable changes in bubble and dew point pressures, as well as critical condensation pressures and temperatures. CO2 injection further complicates the system, enhancing interactions and expanding the coexistence region of the liquid and gas phases on the pressure–temperature diagram, especially across varying pore sizes. Optimization research on CO2 huff and puff technical parameters for shale oil reservoirs suggests the following optimal settings: a CO2 injection rate of 100 t/day, a shut-in time of 30 days, and six huff and puff cycles. The results of this study offer critical insights into CO2-EOR mechanisms in shale oil reservoirs and emphasize the importance of nanopore properties in EOR.

A dynamic Pre-Darcy model and its application in the numerical simulation of enhanced geothermal system

Numerical simulations of Enhanced Geothermal Systems (EGS) currently primarily depend on the Darcy model, neglecting pre-Darcy flow in the matrix rock and the impact of temperature on pre-Darcy flow. This paper introduces a dynamic pre-Darcy model that accounts for the influence of temperature, deriving from the capillary flow model and boundary layer theory. This model is able to consider how temperature variations affect both the thickness of the boundary layer and the yield stress of the fluid. The reliability of the dynamic pre-Darcy model is validated through experimental data, and comparison with both static pre-Darcy and Darcy models demonstrates its superior ability to accurately capture the influence of temperature on pre-Darcy flow. Subsequently, a thermal-hydraulic simulation methodology incorporating the dynamic pre-Darcy model is developed to characterize mass and heat transfer between the matrix and fracture based on the framework of the embedded discrete fracture model (EDFM), which is further validated with analytical and numerical solutions. In the simulation of single-fracture case, we explore the impact of the dynamic pre-Darcy flow on mass and heat transfer between the matrix and fracture. Field-scale numerical simulations for EGS further demonstrate that the Darcy model tends to overestimate both production temperature and total extracted energy. In contrast, the static pre-Darcy model generally underestimates these metrics. The dynamic pre-Darcy model, however, provides more accurate and reasonable predictions.

A Novel Connection Element Method for Multiscale Numerical Simulation of Two-Phase Flow in Fractured Reservoirs

Summary This paper presents a novel approach to the numerical simulation of fractured reservoirs, called the connection element method (CEM), which differs from traditional grid-based methods. The reservoir computational domain is discretized into a series of nodes, and a system of connection elements is constructed based on the given connection lengths and angles. The pressure diffusion term is approximated using generalized finite difference theory. Meanwhile, the transmissibility and volume of the connection elements are determined, and pressure equations are solved discretely to obtain pressure at nodes to approximate the upstream flux along connection elements. Then, we solve the transport equation to obtain oil saturation profiles with low numerical diffusion, utilizing the discontinuous Galerkin (DG) method. Moreover, the flow path tracking algorithm is introduced to quantify the flow allocation factors between wells. In all, the pressure equation can be solved at a global coarse-scale point cloud and the saturation equation is calculated at a local fine-scale connection element. In other words, CEM is of multiscale characteristics relatively. Finally, several numerical examples are implemented to demonstrate that CEM can achieve a relatively better balance between computational accuracy and efficiency compared with embedded discrete fracture modeling (EDFM). Furthermore, CEM adopts flexible meshless nodes instead of grids with strong topology, making it more practical to handle complex reservoir geometry such as fractured reservoirs.

A five field formulation for flow simulations in porous media with fractures and barriers via an optimization based domain decomposition method

The present work deals with the numerical resolution of coupled 3D–2D problems arising from the simulation of fluid flow in fractured porous media modeled via the Discrete Fracture and Matrix (DFM) model. According to the DFM model, fractures are represented as planar interfaces immersed in a 3D porous matrix and can behave as preferential flow paths, in the case of conductive fractures, or can actually be a barrier for the flow, when, instead, the permeability in the normal-to-fracture direction is small compared to the permeability of the matrix. Consequently, the pressure solution in a DFM can be discontinuous across a barrier, as a result of the geometrical dimensional reduction operated on the fracture. The present work is aimed at developing a numerical scheme suitable for the simulation of the flow in a DFM with fractures and barriers, using a mesh for the 3D matrix non conforming to the fractures and that is ready for domain decomposition. This is achieved starting from a PDE-constrained optimization method, currently available in literature only for conductive fractures in a DFM. First, a novel formulation of the optimization problem is defined to account for non permeable fractures. These are described by a filtration-like coupling at the interface with the surrounding porous matrix. Also the extended finite element method with discontinuous enrichment functions is used to reproduce the pressure solution in the matrix around a barrier. The method is presented here in its simplest form, for clarity of exposition, i.e. considering the case of a single fracture in a 3D domain, also providing a proof of the well posedness of the resulting discrete problem. Four validation examples are proposed to show the viability and the effectiveness of the method.

Numerical analysis of coupled thermal–hydraulic processes based on the embedded discrete fracture modeling method

The enhanced geothermal system (EGS) technology is considered an effective way to extract energy from deep hot dry rock (HDR) geothermal resources, which involves the coupling of complex thermal–hydraulic processes in fractured geological media. In this work, we propose a numerical approach based on the concept of embedded discrete fracture modeling method (EDFM), to estimate heat production performance from the EGS. We validated the method by comparing the outputs of our numerical models with experimental data and analytical solutions. A simplified, idealized two-dimensional model with discrete fractures was applied to gain insights on mass and heat transfer mechanisms expected under field conditions. Simulation results indicated that isolated and dead-end fractures have little influence on the overall performance of EGS, under the initial and boundary conditions of our numerical experiments (e.g., fracture density and permeability differences between fractures and matrix). It emerges that geometry and hydraulic transmissivity of fractures are key factors in controlling the effectiveness of the heat exchange between host rocks and circulating fluids. This work emphasizes the importance of implementing increasingly reliable descriptions of the fluid-dynamic and thermodynamic behavior of fractures, to increase the confidence in the numerical assessment of the performance of EGS installations.

Effect of Cross-Well Natural Fractures and Fracture Network on Production History Match and Well Location Optimization in an Ultra-Deep Gas Reservoir

Understanding subsurface natural fracture systems is crucial to characterize well production dynamics and long-term productivity potential. In addition, the placement of future wells can benefit from in-depth fracture network connectivity investigations, vastly improving new wells’ profitability and life cycles if they are placed in dense, well-connected natural fracture zones. In this study, a novel natural fracture calibration workflow is proposed. This workflow starts with the extraction of sector geology and a natural fracture model from the pre-built full-field model. Then, a cross wellbore discrete fracture network (CW-DFN) is created using a novel CW-DFN generation tool, based on image log data. An innovative fracture network identification tool is developed to detect the interconnected regional fracture network (IcFN) with CW-DFN. The non-intrusive embedded discrete fracture model (EDFM) is utilized to numerically incorporate the complex IcFN and CW-DFN in a reservoir simulation, and it is history-matched by tuning their conductivities. This workflow is applied to a single vertical well within a natural fracture carbonate reservoir in Northwest China. The study results show that the number of CW-DFNs is 11, and the number of IcFNs is 72. The non-intersected natural fractures only account for 5.5% of the production, and thus can be removed to improve simulation efficiency. The history-matching absolute average relative deviation (AARD) is 15.16%. The calibrated effective fracture permeability is 280 millidarcy, with an aperture of 0.001 m, equating to a conductivity of 0.28 millidarcy-meter. The 30-year gas production forecast is estimated to be 1.66 billion cubic meters based on a history-matched model. Finally, if the well is drilled to the east of the sector, 30-year production declines to 1.33 billion cubic meters (a reduction of 20%). However, if the well is drilled to the west of the sector, 30-year production increases to 2 billion cubic meters (an improvement of 20.5%).

Coal particle size-dependent gas diffusion coefficient

The tests on gas (CH4 or CO2) desorption-diffusion in coal particles, together with the analytical unipore or bidisperse diffusion model, are widely used to study the gas diffusion behavior, but the fitted diffusion coefficient always increases with particle size, which makes it difficult to determine an identical diffusion coefficient for coal particles. In this work, five sets of anthracite particles with different sizes are prepared, and tests on CO2 desorption-diffusion in different coal particles are carried out. Four continuum models (including an analytical unipore diffusion model, an analytical bidisperse diffusion model, a numerical unipore diffusion model, and a dual-porosity model) and a discrete fracture model are established. The continuum models can all fit the experimental results well, while the fitted diffusion coefficients always increase with particle size. Based on the identical matrix’s side length, matrix diffusion coefficient, and fracture diffusion coefficient, the discrete fracture model fits the experimental results of different coal particles successfully, and the size effect of diffusion coefficient is resolved. The diffusion coefficient fitted by the continuum models is the equivalent diffusion coefficient of matrix diffusion and fracture diffusion, so the equivalent diffusion coefficient increases with fracture path (i.e., particle size), and the size effect of diffusion coefficient occurs. The results provide an insight into the size effect of diffusion coefficient and make it possible to determine an identical diffusion coefficient for coal particles of different sizes.

Physics-informed neural network simulation of two-phase flow in heterogeneous and fractured porous media

Physics-informed neural networks (PINNs) have received great attention as a promising paradigm for forward, inverse, and surrogate modeling of various physical processes with limited or no labeled data. However, PINNs are rarely used to predict two-phase flow in heterogeneous and fractured porous media, which is critical to lots of subsurface applications, due to the significant challenges in their training. In this work, we present an Enriched Physics-Informed Neural Network (E-PINN) to overcome these barriers and realize the simulation of such flow. Specifically, the Embedded Discrete Fracture Model (EDFM) is adopted to explicitly represent fractures, and then the finite volume method (FVM) instead of the Automatic Differentiation (AD) is used to evaluate spatial derivatives and construct the physics-informed loss function, so that the flux continuity between neighboring elements with different properties (e.g. matrix and fracture) can be defined rigorously. Besides, we develop a novel physics-informed neural network (NN) architecture adopting the adjacency-location anchoring, adaptive activation function, skip connection and gated updating to enrich the pressure information and enhance the learning ability of NN. Additionally, the initial and boundary conditions are constrained through a hard approach, which encodes them into network design, to improve the accuracy and efficiency of network training. In order to further reduce the difficulty of training, the Implicit-Pressure Explicit-Saturation (IMPES) scheme is used to calculate pressure and saturation, in which only the pressure needs to be solved by training NN. Finally, the superiority and applicability of E-PINN to complex practical problems is demonstrated through the simulations of immiscible displacement in 2D/3D heterogeneous and fractured reservoirs.

Complex fracture description and quasi-elastic energy development mathematical model for shale oil reservoirs

Reservoir numerical simulation is an important tool and method for the reasonable and efficient development of shale reservoirs. Accurate description of three-dimensional fractures in shale reservoir development is a necessary and sufficient condition to improve the accuracy and robustness of shale reservoir numerical simulation. This paper achieves precise characterization of complex fracture shapes and oil, gas and water flow by establishing an embedded discrete fracture model based on a non-structural network, which has advantages in the fine characterization of complex morphological fractures in the reservoir and the grid division of the reservoir. In the large matrix solution method, the Newton-Raphson method is used to linearize the nonlinear equations, the Jacobian matrix is ​​constructed, the ILU method is used for preprocessing, the conjugate gradient method is used to solve the linear equations, and the shale oil quasi-elasticity is established A fully implicit solution method for mathematical models of energy development.

An extended isogeometric collocation method for fracture analysis

AbstractA collocation method is developed for discrete fracture models in the context of the partition‐of‐unity method. Spline technologies used in isogeometric analysis (IGA) are exploited to provide a smooth inter‐element transition of gradients, thus allowing to get rid of extra flux terms at element boundaries which are generated by Lagrange polynomials. Bézier extraction is utilised to formulate IGA commensurate with a standard finite element data‐structure. The efficacy of the proposed approach is examined through different numerical examples and is compared with other discrete methods for fracture analysis. The proposed approach is competitive in terms of accuracy with the least computational cost, rendering it a suitable candidate for superseding available collocation approaches for fracture simulation. Moreover, the approach naturally assesses the possibility of physics informed neural networks for fracture simulation, to which collocation is central.

A new analytically modified embedded discrete fracture model for pressure transient analysis in fluid flow

In the past few decades, multi-stage hydraulic fracturing has emerged as a crucial technology for the commercial development of unconventional oil and gas (UOG) resources. It is crucial to accurately and efficiently characterize transient fluid flow near fractures, which is a critical concern for many researchers. Discrete fracture models (DFMs) are primarily used to analyze the pressure transient behaviors of fluid flow in naturally fractured porous media. DFMs can accurately capture transient fluid flow near fractures, but they require a substantial number of grids to ensure computational accuracy, leading to higher computational costs. On the other hand, standard embedded discrete fracture models (EDFMs) based on pseudo-steady-state assumptions are computationally efficient, but they struggle to model the early transient fluid flow near fractures accurately. To address this limitation, we propose a new analytically modified EDFM (AEDFM) with structured Cartesian grids for analyzing the pressure transient behaviors of fluid flow in naturally fractured porous media.The transmissibility between the matrix and fractures is adjusted by multiplying it with a transient factor. In addition, we have validated the accuracy and efficiency of our proposed model through comparisons with results from analytical models and a standard well-testing software package. The results demonstrate the significance of our proposed model in accurately capturing transient fluid flow around fractures and reducing computational costs. In this work, we analyze the pressure transient behaviors of fluid flow using various parameter values and further evaluate the significance of the proposed modifications. The results indicate that AEDFM can effectively match the early nonlinear pressure drop near fractures compared to the standard EDFM. This work presents a powerful tool for the fast and accurate analysis of pressure transient behaviors of fluid flow in naturally fractured porous media.

Influence of fractures and heat transmission on wormhole propagation in fractured carbonate rocks

Carbonate reservoirs exhibit significant heterogeneity and often contain natural fractures. Matrix acidizing is a commonly used technique in fractured carbonate formations. Currently, the influences of natural fractures and the acid-rock molar reaction heat on wormhole propagation are rarely considered in acidizing models. In this paper, the two-scale continuum model, the embedded discrete fracture models (EDFM), and the heat transfer model are combined to develop a computational model of acidic fluid reactive flow in carbonate rocks with a complex fracture network that takes into account temperature variations. The model considers not only the temperature variation caused by the acid-rock molar reaction heat, but also simulates carbonate acidizing with any complex fracture network. Numerical simulations of the model were carried out to obtain the wormhole propagation patterns and the breakthrough curves under different conditions, and the simulation results were in good agreement with those of previous study. In particular, the sensitivity of fracture permeability, acid concentration, acid temperature and matrix porosity heterogeneity during acidizing is studied. The simulation results show that fractures can significantly alter the structure of wormholes and decrease the pore volume of acid required to reach the breakthrough (PVBT). However, fractures have minimal impact on the growth structure of wormholes prior to encountering fractures. During the reaction, the highest temperature was observed at the front end of the wormhole, and it was found that the optimum injection rate increased with the increase of the injected acid temperature. Lowering the temperature of the acid can enhance the efficiency of acidizing in fractured carbonate formations, particularly when applied at an optimal injection rate. When carbonate formations have well-developed fractures, the structure of the wormholes is insensitive to the heterogeneity of the matrix porosity, and they tend to grow along the fractures.

Meshfree multiscale method for the infiltration problem in permafrost

This study introduces a multiscale simulation approach, employing the Meshfree Generalized Multiscale Finite Element Method (GMsFEM), to numerically model fluid seepage in permafrost soil conditions. Combining Meshfree and Finite Element Methods, the coarse-grid system is addressed, while the standard Finite Element Method operates on the fine grid. For representing fractures at the fine grid level, the Discrete Fracture Model (DFM) is utilized. The findings demonstrate the method’s efficacy in accurately depicting soil heterogeneity and the impact of fractures on soil temperature.

A Practical probabilistic history matching framework for characterizing fracture network parameters of shale gas reservoirs

Simulating fluid flow in complex fracture systems encountered in subsurface systems, such as those encountered in tight or shale gas reservoirs, presents considerable challenges. Sophisticated numerical techniques and characterization approaches are required to integrate static and dynamic data from multiple sources to account for the intricate interplay between fractures and the surrounding rock matrix. A novel workflow using an indicator-based probability perturbation method (PPM) is applied to characterize uncertain discrete fracture network parameters. The posterior probability distributions of primary and secondary fracture transmissivity, aperture, length, height, and intensity of the secondary fracture are estimated according to the production (flow) histories. A case study based on the Horn River shale gas reservoir is presented. A pilot point scheme is formulated to update the distributions of P32L. The forward modelling entails upscaling each fracture model into a dual-porosity dual-permeability model and performing multiphase flow simulation. This approach produces an ensemble of history-matched discrete fracture and upscaled models by incorporating static and dynamic data. The results demonstrate the utility of the developed method for estimating secondary fracture parameters, which are not inferrable from other static information alone. The inference of both primary and secondary fractures offers a more comprehensive characterization of the fracture networks in tight or shale gas reservoirs. Integrating a pilot point parameterization technique and a sequential simulation approach into the PPM framework presents significant novelties in contributing to a comprehensive and effective method of simulating fluid flow in complex fracture systems, addressing the challenges associated with non-linear relationships between fracture model parameters and the corresponding flow response.

A novel multiphase and multicomponent model for simulating molecular diffusion in shale oil reservoirs with complex fracture networks

Molecular diffusion is critical for enhanced oil recovery (EOR) in shale oil reservoirs with complex fracture networks. Understanding the influence of fractures on diffusive mass transfer is crucial for predicting oil recovery and remaining oil distribution. Diffusive mass transfer between fractures and matrix is critical in comprehensively and effectively simulating molecular diffusion. Resolution of matrix cells significantly affects diffusion accuracy at the fracture–matrix interface. Low resolution results in multiple fractures in the same matrix cell, leading to decreased precision in calculating mass transfer by conventional methods. To address this, a novel multiphase and multicomponent model is proposed. The new model integrating the consideration of fracture spacing modifies molecular diffusion transmissibility between fracture and matrix in an embedded discrete fracture model. The discretization employs the two-point flux approximation in the finite-volume method. Validation compares the coarser mesh to the finest grid as a reliable reference. Results show the proposed model accurately captures diffusive mass transfer in a coarser mesh. Modified models study molecular diffusion's effects on EOR in shale oil reservoirs with complex fracture networks by CO2 huff and puff. Results indicate that increasing injection rates cannot improve oil recovery under extremely low porosity and permeability. Molecular diffusion facilitates CO2 penetration into the formation. This expands the swept CO2 volume and increases both volume expansion and formation energy. In addition, the light and heavy components of the crude oil are diffused into the fractures and eventually produced, which reduces gas production in the case of diffusion.

A New Simplified Discrete Fracture Model for Shearing of Intersecting Fractures and Faults

Shearing of fractures and faults is important because it can result in permeability change or even induce seismicity—both are keys for efficient and safe energy recovery and storage in Earth systems. Quantitative analysis of shearing of intersecting fractures and faults is challenging because it can involve dynamic frictional contacts that are complicated by deformation of the rock matrix. To predict the shearing of intersecting fractures/faults, we attempt to answer the question of how intersections impact the shearing of a fracture network and whether we can simplify the description as compared to classical discrete fracture network (DFN) models. To answer these questions, we conducted a series of numerical simulations on scenarios for variable numbers of intersecting fractures. All these examples yield consistent results: the results of using DFNs are consistent with those of using hypothetical major paths. This leads to a new model, which we name simplified discrete fracture network model, to analyze shearing of intersecting fractures/faults using major path(s). We found that the intersections of fractures do not fundamentally change the shearing of two intersecting fractures if the intersecting angles are small. Furthermore, increasing the number of fractures/faults may relax the stress as more fractures/faults become available for shearing and distributing the stress. The simplified DFN model, which can capture efficiently the shearing behavior of each major paths from a large number of intersecting fractures/faults, will be a promising conceptual model that is complementary to existing equivalent continuum and discrete fracture models to analyze shearing of intersecting fractures/faults.

Numerical methods for reinterpreted discrete fracture models with random inputs

This paper investigates the impact of uncertainty on the reinterpreted discrete fracture model (RDFM) for flow in porous media featuring fractures and barriers. The stochastic RDMF is formulated by parameterizing the uncertainty in the permeability, as well as the width of the fracture and barrier networks as random variables. To discretize these random variable, we adopt a stochastic collocation method based on a special choice of nodes. This collocation scheme offers optimal accuracy and easy implementation. The resulting collocation scheme is then combined with the local discontinuous Galerkin (LDG) method to accommodate the hybrid-dimensional Darcy’s law and the nature of the pressure/flux discontinuity. Additionally, a slope limiter is utilized to enhance stability. Numerical experiments are carried out to illustrate the accuracy and stability of the proposed scheme.