# A New Multicontinuum Model for Advection-Diffusion Process of Single-Phase Nonlinear Flow in a Multiscale Fractured Porous Media

Authors: Richard Owusu; Adu Sakyi; Peter Amoako-Yirenkyi; Isaac Kwame Dontwi
Journal: Journal of Applied Mathematics
 (2022)

---

## Abstract

Fractured porous media modeling and simulation has seen significant development in the past decade but still pose a great challenge and difficulty due to the multiscale nature of fractures, domain heterogeneity, and the nonlinear flow fields due to the high flow velocity and permeability resulting from the presence of fractures. Therefore, modeling fluid transport that is influenced by both advection and diffusion in fractured porous media studies becomes a generic problem, which this study seeks to address. In this paper, we present a study on non-Darcian fluid transport in multiscale naturally fractured reservoirs via an upscaling technique. An averaged macroscopic equation representing pressure distribution in a three-phase multiscale fractured porous medium was developed, consisting of the matrix and a 2-scale fractured network of length-scalesℓm and ℓM. The resulting macroscopic model has cross-advective and diffusive terms that account for induced fluxes between the interacting domains, as well as a mass transfer function that is dependent on both physical and geometric properties of the domain, with both advective and diffusive properties. This model also has effective diffusive and advective coefficients that account for reservoir properties such as viscosity, fluid density, and flow velocity. From the numerical simulation, a radial, a horizontal-linear flow behavior, and a transient and quasi-steady-state flow regime that is typical of naturally fractured porous media was observed. The findings of this study will provide researchers a reliable tool to study fractured porous media and can also help for better understanding of the dynamics of flow in fractured reservoirs.


---

## Body

## 1. Introduction

Natural porous media such as fractured reservoirs are not uniformly random or homogeneous in general. It may contain multiple fractures or fissures of varying sizes. In petroleum reservoirs, the presence of cracks, fissures, and fractures can cause a high degree of heterogeneity in the porous media, resulting in different flow properties occurring at different points in the domain of varying fracture lengths or sizes. According to Altinörs and Önder (2010) [1], a nonhomogeneous porous medium creates a nonuniform flow field that affects movement of fluids with significant different flow velocities and flow behavior [1, 2], a common phenomena identified with multiscale fractured porous media. The effective representation of advection and diffusion in models associated with fluid transport in naturally fractured reservoirs is a generic problem for researchers and oil field engineers. Because of the multiscale heterogeneity in the underlying mediums, modeling advection and diffusion transport in these flow fields become a complex process.To successfully model these flow processes, a more tactful approach is required. Various attempts have been made via the multicontinuum approach, resulting in triple-porosity and triple-porosity/permeability models ([1, 3–7]) that are superior to dual porosity and permeability models ([3, 8]). However, because they are based on much simplified assumptions, these models are less reliable to study advection/diffusion in fractured porous medium with nonlinear flow behaviors. Among the limitations of these models are laminar flow assumptions and a less accurate representation of mass transfer functions for fracture modeling [4]. Furthermore, representing macroscopic flow equations with original mass conservation equations and Darcy equations results in models that lack important flow properties such as cross advection/diffusion effects, with inefficient representation of diffusivity coefficients and the mass transfer function. Cross advection and diffusion have been identified and extensively studied with significant effects in reactive diffusive systems ([9, 10]), in which diffusion and advection in one domain are induced by other interacting domains, typical of naturally fractured porous medium.To obtain an efficient transport model for advection and diffusion, macroscopic flow equations derived from point microscopic transport equations via an upscaling technique are recommended ([4, 11]). Volume averaging method (VAM) is a homogenization or an upscaling technique that derives macroscale variables through spatial averaging of heterogeneous media. Through the application of VAM, macroscale equations with closure variables can be derived. This leads to solving a heterogeneous problem at the macroscale level that represents the original heterogeneous medium. In VAM, averaged quantities of the original variables are used to preserve the original properties of the variables. The main objective of this study is to develop a local-volume averaged transport model with effective flow properties that depend on some lower-scale complex heterogeneous medium for fluid transport in a multiscale naturally fractured porous medium.This goal is primarily achievable through the use of upscaling techniques, which are capable of representing microscopic flow properties at the macroscopic level. The method allows us to model fluid transport in the proposed heterogeneous fractured porous medium, consisting of the matrix and the 2-scale fractured porous network with fracture length-scalesℓm and ℓM as shown in Figure 1. Using the volume averaging concept, we were able to model flow interactions in the matrix as well as each of the microscale and macroscale fracture networks. As a result, three local averaged transport equations were developed, one for fluid transport in the matrix domain and the other two for fluid transport in the 2-scale fractured network domain.Figure 1
Proposed multiscale naturally fractured porous media model.From our findings, our proposed model yielded a mass transfer function having both convective and diffusive property and additional microscopic reservoir properties including flow velocity, fluid viscosity, compressibility, and surface normal vector for fracture orientations, which are lacking in the commonly used transfer function originally proposed by Barenblatt et al. [3] and Warren and Root [2]. The upscaling technique employed in the study resulted into additional microscopic properties such as cross-advection and cross-diffusion terms, which account for induced fluxes due to flow interaction between the matrix and the fractures. Comparatively, existing models [8, 12] fail to account for this unique features obtained in our model, which makes our proposed model more reliable for porous media study. In addition, this model incorporates nonlinear flow behavior which enhance the suitability of the model for fractured porous media study, a characteristic which is identical to fractured porous media, but lacking in existing models.The study began with the fundamental mass conservation and Darcy-Forchheimer momentum equations. A nonlinear transport model obtained from the fundamental equations was used to set-up the model problem at the microscopic level, considering the possible interactions at the matrix-fracture interfaces. Then, the microscopic model was volume averaged to obtain a volume averaged transport equation (VATE), from which a quasi-steady closure problem was developed for some spatial deviation variables obtained in the VATE. The closure problem was then solved for the deviation variables, which were used to obtain a closed form problem of the VATE model at the macroscopic level.Section2 gives the governing equations and the theories of volume averaging. In Section 3, we exhibit the spatial smoothing approach for two-phase heterogeneous media where a closure problem is developed and a solution is proposed for the deviation variables. Sections 4 and 5 present the model proposed in the study. Finally, we discussed the results of the study in Section 6 and concluded in Section 7.


## 2. Governing Equations and Volume Averaging Theorems

The fractured porous medium considered in this study is assumed to consist of the matrix and a 2-scale fractured porous network with fracture length-scalesℓm and ℓM with different porosity and permeability, hence different flow characteristics. The set up of the model (1) considers mass transfer between the matrix and each of the two scales of fractures and (2) assumes flow to be nonlinear due to the high flow velocity and permeability that exist in a fractured porous medium ([1]). The governing equations of motion describing fluid transport in the fractured porous domain at the microscopic level are represented by the mass conservation
(1)∂tΦρ+∇·ρq=0inD×J,where Φ, ρ, and q represent porosity, fluid density, and flow velocity, and the Forchheimer equation
(2)−∇P=μKq+βρqq,which relates the nonlinear pressure drop to the flow velocity for porous media with high flow velocity and permeability, where P, μ, K, and β, respectively, denote pressure, fluid viscosity, domain permeability, and Forchheimer coefficient that depends on permeability and porosity. Secondly, at the microscopic scale level, incompressible flow is assumed thus, ∇·q=0.Consider an averaging volume with two phasesα and β, as presented in Figure 2. Given the quantity Ψ in phase α, the volume average for an averaging volume of radius r is defined by.
(3)Ψα=1V∫ε∈VxmεΨαx,tdVε,or using the intrinsic properties of the averaging volume, the intrinsic volume average is obtained as follows:
(4)Ψαα=1Vα∫ε∈VαxmεΨαx,tdVε,where
(5)mε=1ifε≤r0ifε>rand Ψα and Ψαα are related by
(6)Ψα=εαΨαα,εα=VαV.Figure 2
A two possible averaging volumes [13].Volume averaging the spatial derivative of the quantity in phaseα is obtained by
(7)∇Ψα=1V∫Sαexnα·ΨαdS+1V∫Sαinnα·ΨαdS,which is simplified as follows [13]:
(8)∇Ψα=∇Ψα+1V∫Sαinnα·ΨαdS.


## 3. Spatial Smoothing for Two-Phase Heterogeneous Media

Consider a two-phase heterogeneous system, consisting of anα- and β-phase. For each point in the domain, we define a microscale averaging volume having internal heterogeneities with permeable interfaces. The governing equations for fluid transport were developed from the mass conservation and momentum equations defined in equations (1) and (2) [7]. Thus, an initial-boundary value problem for fluid flow in α-phase in an averaging volume, V, is obtained as follows:
(9a)cαΦα∂tPα+∇·cαqαPα−∇·Kα∇Pα=0inphaseα,(9b)∇·qα=0inphaseα,(9c)BC1:−Kα∇Pα·nαβ=kαPαatAαβ,(9d)BC2:Pαx,t=fx,tatAαe,(9e)IC:Pαx,0=fxatt=0,where cα denotes compressibility, Kα is the diffusive coefficient in phase α, kα represents the pseudo reaction rate constant for heterogeneous reaction at the α-β interface, Aαβ is the area of the α-β interface contained within the averaging volume, Aαe denotes the exits and entrances of the α-phase at the boundaries of the macroscopic region, and V is expressed as follows:
(10)V=Vα+Vβ,where Vα is the α-phase volume contained within the averaging volume V expressed as εα=Vα/V and the BC1 is imposed to account for interface reaction and flux contribution at Aαβ [14]. To carry out the spatial smoothing process, we define the superficial and intrinsic average quantities
(11)Pα=1V∫VαPαdV,Pαα=1Vα∫VαPαdV,of the pressure variable Pα, respectively, which are related by
(12)Pα=εαPαα.To begin with, the point advective-diffusion equation (9a) is spatially smoothed and the fact that the averaging volume is time invariant leads to
(13)cαΦα∂tPα+∇·cαqαPα−∇·Kα∇Pα=0.Using the spatial average theorem for the second and third terms in equation (13) leads to the averaged transport equation simplified as
(14)Φαcα∂tPα+∇·cαqαPα−∇·Kα∇Pα+1V∫AαβnαβPαdA=−1V∫AαβkαPαdA−cαV∫Aαβnαβ·qαPαdA,where the interfacial BC1 in (9c) is used and the integral term
(15)1V∫AαβkαPαdA=SvkαPα,Sv=AαβV,Pα=1Aαβ∫AαβPαdA,is referred to as the heterogeneous reaction term accounting for internal heterogeneities [14, 15]. Since we are estimating the quantity, Pα, by an average quantity Pα, we account for deviations by incorporating the spatial deviation variable Pα~ according to the relation:
(16)Pα=Pα+Pα~.Substituting identity (16) into (13) gives
(17)cαΦαεα∂tPαα+∇·cαqαPα−∇·Kαεα∇Pαα+1V∫AαβnαβPα~dA=−SvεαkαPα−cαV∫Aαβnαβ·qαPα~dA,where the terms with the volume averaged quantity, Pα in the area integrals
(18)1V∫AαβnαβPαdA,cαV∫Aαβnαβ·qαPα,are removed (nonlocal terms) [14, 15] and the intrinsic average for Pα represented. To analyze the convective term, ∇·cαqαPα, the spatial deviations of Pα and qα are expressed as follows:
(19)Pα=Pαα+Pα~,qα=qαα+qα~,from which qαPα can be expanded as
(20)qαPα=qαα+qα~Pαα+Pα~=qααPαα+Pααqα~+Pα~qαα+qα~Pα.~By neglecting variation of averaged quantities within the averaging volume (i.e.,Pαα=Pαα) and since the average of the spatial deviation Pα~α=0 [16], we have that
(21)cαqαPα=cαεαqααPαα+cαqα~Pα~.Substituting (21) into (17) results in the volume averaged transport equation (VATE)(22)cαΦαεα∂tPαα+∇·cαεαqααPαα−∇·Kαεα∇Pαα+1V∫AαβnαβPα~dA=∇·cαqα~Pα~−SvεαkαPαα−cαV∫Aαβnαβ·qαPα~dA,where the terms ∇·cαεαqααPαα and ∇·cαqα~Pα~ are the advective transport and dispersive transport, respectively. To obtain a close form problem at the macroscopic level, explicit representation of Pα~ is found in terms of the averaged quantity Pαα by developing a boundary value problem for the spatial variables. The process begins by dividing through equation (22) by εα and subtracting from the original
(23)cαΦα∂tPα+cαqα·∇Pα−Kα∇·∇Pα=0,which is further simplified to obtain a partial differential equation for the spatial deviation variable as follows:
(24)cαΦα∂tPα~⏟Accumulation+cαqα·∇Pα~⏟advection+cqα~·∇Pαα⏟advectionsource=∇·Kα∇Pα~⏟diffusion+εα−1∇εα·Kα∇Pαα⏟+εα−1∇·KαV∫AαβnαβPα~dA⏟non−localdiffusion+εα−1∇·cαqα~Pα~⏟non−localadvection−SvkPαα−εα−1cαV∫Aαβnαβ·qαPα~dA⏟reactivesourcediffusionsource.The differential equation (24) is to be solved at the microscopic level and not the macroscopic level; hence, noncontributing terms at the microscopic level are identified and removed. To do this, we let L and ℓα be the characteristic length scale of the macroscopic medium and the length scale of the α-phase of the averaging volume which satisfies constraint
(25)ℓα<<L.We estimate the spatial deviation variablePα~ in the order of the averaged quantity Pα as follows:
(26)Pα~=OPα,and thus
(27)∇Pα=O∇PαL,∇Pα~=O∇Pαℓα.From these estimates, the nonlocal diffusion and the diffusion source terms are discarded in comparison to the diffusion term while the nonlocal advection term is discarded when compared to the advection term. The mass accumulation term,cαΦα∂tPα~, is discarded under quasi-steady conditions, which finally leads to the quasi-steady closure problem expressed as follows [14, 17]:
(28)−∇·Kα∇Pα~+cαq·∇Pα~=−cαqα~·∇Pαα−SvkαPαα−εα−1cαV∫Aαβnαβ·qαPα~dABC1:−nαβ·Kα∇Pα~=nαβ·Kα∇Pα+kαPαatAαβBC2:Pα~x+li=Pα~xi=1,2,3,where the BC2 in equation (28) is replaced by a local condition (periodic condition), since we are to solve the closure problem in a unit cell representation that is spatially periodic at the microscopic level.It is obvious that the spatial deviation variable will be dependent on∇Pαα and Pαα. Thus, Pα~ can be expressed in terms of the closure variables b and s as [13, 14, 18]
(29)Pα~=b·∇Pαα+sPαα,which is used together with equation (22) to derive a closed form averaged transport equation at the macroscopic level.


## 4. Advection and Diffusion for Non-Darcian Flow in a Three-Phase Multiscale Fractured Porous Media System

Consider a single-phase fluid transport in naturally fractured porous domainD, comprising of the α-, β-, and γ-phase, which represent the matrix and the micro- and macrofracture networks, respectively (see Figure 3). Here, we consider a single-phase fluid transport process in the α,β,γ-system illustrated in Figure 3. The governing equations for the proposed multiscale heterogeneous fractured porous medium are given as follows (in the α-region):
(30)cαΦα∂tPα+cα∇·qαPα−∇·Kα∇Pα=0inphaseα,BC1:Pαx,t=Pηx,tatAαη,BC2:Kα∇Pα·nαη=Kη∇Pη·nαηatAαη,and in the η-region (η∈β,γ) as follows:
(31)cηΦη∂tPη+cη∇·qηPη−∇·Kη∇Pη=0inphaseη,BC3:Pαx,t=fx,tatAαe,BC4:Pηx,t=gx,tatAηe.Figure 3
Fractured reservoir model: a macroscopic region and an averaging volume [7].For a three-phase system considered in this study, we express the averaging volume,V, shown in Figure 3 as follows:
(32)V=Vα+Vβ+Vγ,where the volume fraction, εi for phase i∈α,β,γ, in the averaging volume is expressed as
(33)εi=ViVi∈α,β,γ.Following similar analysis for spatial smoothing in Section3, we obtain the volume averaged transport equation for the α-η system written as follows (in α-region):
(34)cαΦαεα∂tPαα+∇·cαεαqααPαα−∇·Kαεα∇Pαα+1V∫AαηnαβPα~dA=∇·cαqα~Pα~+1V∫Aαηnαβ·Kα∇PαdA−cαV∫Aαηnαβ·qαPα~dAand in the η-region we have
(35)cαΦηεη∂tPηη+∇·cηεηqηηPηη−∇·Kηεη∇Pηη+1V∫AηαnηαPη~dA=∇·cηqη~Pη~+1V∫Aηαnηα·Kη∇PβdA−cηV∫Aηαnηα·qηPη~dA,where the interfacial flux contributing to mass transfer from α-region into η-region is simplified as [19]
(36)1V∫Aαηnαη·Kα∇PαdA=SαηKαPαα−Pηη,and Sαη is the total interfacial surface area at the α,η-interface. The closure problem is obtained for Pα~ and Pη~ following similar analytic technique in Section 3 and written as in the α-region
(37)−∇·Kα∇Pα~+cαqα·∇Pα~=−cαqα~·∇Pαα+Sαηεα−1KαPαα−Pηη−εα−1cαV∫Aαηnαη·qαPα~dA,BC1:Pα~=Pη+Pη~−PαatAαη,BC2:nαη·Kα∇Pα~=nαη·Kα∇Pη+nαη·K∇Pη~−nαη·Kα∇PαatAαη,and in the η-region, we have
(38)−∇·Kη∇Pη~+cηqη·∇Pη~=−cηqη~·∇Pηη+Sηαεη−1KηPηη−Pαα−εη−1cηV∫Aηαnηα·qηPη~dA,BC3:Pα~x,t=fx,tatAαe,BC4:Pη~x,t=gx,tatAηe.Considering the closure problem, the source terms and the boundary terms are represented in terms of the averaged variables identified as∇Pαα,∇Pηη,andPβ−Pα. It is therefore obvious for the solution of the two-equation model to be dependent on the source terms. Hence, following the concept in the works of [7, 17, 20, 21], the spatial deviations are represented as functions of the averaged variables of the macroscopic source terms written as follows:
(39)Pα~=bα·∇Pαα+bαη·∇Pηη+sαPηη−Pαα,Pη~=bη·∇Pηη+bηα·∇Pαα+sηPαα−Pηη.


## 5. The Closure Variables

The closure variables are obtained by substituting the solution proposed in (39) into the closure problem (37) and (38) and simplified, where “non-contributing” terms are removed. This leads to(40)−Kα∇2bα∇Pαα−Kα∇2bαη∇Pηη−Kα∇2sαPηη−Pαα+cαqα∇bα∇Pαα+cαqα∇bαη∇Pηα+cαqαsαPηη−Pαα+Aαbα∇Pαα+Aαbαη∇Pηη+AαPηη−Pαα=−cαqα∇Pαα−SαεαKαPηη−Pααwhere
(41)Aα=εα−1cαV∫Aαηnαηqαη.After this, the right-hand-side and left-hand-side terms in equation (40) are compared, leading to the boundary value problem for the closure variables asProblem 1.
(42)−κα∇2bα+cαqα∇bα+εα−1cαV∫Aαηnαη·qαbα=−cαqαinα−phasebα=bηαnηαKα∇bα−Kη∇bηα=−nηαKα,−κη∇2bηα+cηqη∇bηα+εη−1cβV∫Aηαnηα·qηbηα=0inη−phasePeriodic:bαx+li=bαx,bηαx+li=bηαx.The spatially periodic boundary condition is used in place of the BC3 and BC4 in the closure problem (38). The periodic boundary condition is used in order to solve the problem at the microscopic scale in a unit cell representation. Similarly, the periodic boundary conditions are used in the second and third boundary value problems for the closure variables asProblem 2.
(43)−Kα∇2bαη+cαqα∇bαη+εα−1cαV∫Aαηnαη·qαbαη=0inα−phasebη=bαηnαηKη∇bη−Kα∇bαη=−nαηKη,−Kη∇2bη+cηqη∇bη+εη−1cηV∫Aηαnηα·qηbη=−cηqηinη−phasePeriodic:bηx+li=bηx,bαηx+li=bαηx.Problem 3.
(44)−∇2sα+cαqα∇sα+εα−1cαV∫Aαηnαη·qαsα=−Sαηεα−1Kαinα−phasesα=−sηnαηκη∇sα+Kη∇sη=0,−∇2sη+cηqη∇sη+εη−1cηV∫Aηαnηα·qηsη=−Sηαεη−1Kηinη−phasePeriodic:sαx+li=sαx,sηx+li=sηx.


### 5.1. Closed Form Problem: The Macroscopic Volume Averaged Transport Model

The closed form problem is obtained by substituting equation (39) into equations (34) and (35) and further simplified. This leads to the volume averaged transport equation
(45)cαΦαεα∂tPαα+∇·εαcαqααPαα−∇·εαKα·∇Pαα+Kαη·∇Pβη−∇·εαuαdPβη−Pαα=∇·εαDα·∇Pαα+∇·εαDαη·∇Pηη+∇·uαcPηη−Pαα−SαηKαPηη−Pαα−εαuαs∇Pαα−εαuαβs∇Pηη−εαAvαηPηη−Pαα,where
(46)Kα=KαI+εα−1V∫AαηnαβbαdA,Kαη=εα−1KαV∫AαηnαηbαηdA,uαd=εα−1KαV∫AαηnαηsαdA,Dα=εα−1cαqα~bα,Dαη=εα−1cαqα~bαη,uαc=εα−1cαqα~sα,uαs=εα−1cαV∫Aαηnαη·bαqαdA,uαηs=εα−1cαV∫Aαηnαη·bαηqαdA,Avαη=εα−1cαV∫Aαηnαη·qαsαdA.Equation (45) is further simplified to obtain the macroscopic volume averaged transport equation for fluid flow in the α-region as follows:
(47)cαΦα∂tPαα+∇·qαePαα+uαηsPβη−∇·Kαe·∇Pαα+Kαηe·∇Pβη=∇·uαePβη−Pαα−AvαηPηη−Pαα−SαηKαPηη−Pαα,where
(48)Kαe=Kα+Dα,Kαβe=Kαβ+Dαβ,qαe=cαqαα+uαs,uαe=uαd+uαc,and the right-hand terms ∇·uαePηη−Pαα,AvαηPηη−Pαα and Sαηκα−1Pηη−Pαα represent the mass flux due to advective transport induced by α-η interaction at the interface, interfacial mass transfer by advection, and interfacial diffusive flux contribution at the interacting interfaces, respectively.Again, substitutingPη~ in equation (39) into equation (35) and simplifying leads to macroscopic volume-averaged transport equation in the η-region as follows:
(49)cηΦη∂tPηη+∇·qηePηη+uηαsPαα−∇·Kηe·∇Pηη+Kηαe·∇Pαα=∇·uηePαα−Pηη−AvηαPαα−Pηη−SηαKηPαα−Pηη,with similar representations of the coefficients defined for the α-phase equation. Due to the flow interaction between the α-phase and both the β- and γ-phases, it is important to account for flow interaction effects including cross diffusion and advection effects induced by both β-phase and γ-phase on the α-phase. The cross terms identified in the new equation (47) include uαηsPηη,Kαηe·∇Pηη, and ∇·uαePηη−Pαα where η∈β,γ and other terms arising from flow interaction between the α and both the β and γ phases. Thus, we obtain a system of volume-averaged transport equation for the α,β,γ-systems given as follows:
(50)cαΦα∂tPαα+∇·qαePαα+Uβγ−∇·Kαe·∇Pαα+Dβγ=∇·uαePββ−Pαα+∇·uαePγγ−Pαα−RαβPββ−Pαα−RαγPγγ−Pαα,(51)cβΦβ∂tPββ+∇·qβePββ+uβαsPαα−∇·Kβe·∇Pββ+Kβαe·∇Pαα=∇·uβePαα−Pββ−RβαPαα−Pββ,(52)cγΦγ∂tPγγ+∇·qγePγγ+uγαsPαα−∇·Kγe·∇Pγγ+Kγαe·∇Pαα=∇·uγePαα−Pγγ−RγαPαα−Pγγ,where
(53)Uβγ=uαβsPββ+uαγsPγγ,Dβγ=Kαβe·∇Pββ+Kαγe·∇Pγγ,Rαβ=−Avαβ−SαβKα,Rβα=−Avβα−SβαKβ,Rαγ=−Avαγ−SαγKα,Rγα=−Avγα−SγαKγ.The parameters obtained in this model and the associated parameter values used for the numerical simulation of the model are outlined in Table1. To accurately represent the permeability of the domain, the fracture and matrix permeabilities were evaluated from the relation proposed by [22], which are defined, respectively, as
(54)Kf=0.05758ϕf×hfeVf3f∈β,γ,K=ϕ42s2,where s is the pore area.Table 1
Table of parameters used in the model and parameter values for simulation [12, 23].
ParameterSymbolValue(s)Time (yrs)T5Fluid densityρ1Fluid viscosityμ0.001Porosity of domainϕϕα=0.40, ϕβ=0.05, ϕ=0.04CompressibilityCCα=0.15, Cβ=0.18, Cγ=0.20Fracture apertureeVeVβ=20, eVγ=55Fracture heighthfhβ=100, hγ=115


## 5.1. Closed Form Problem: The Macroscopic Volume Averaged Transport Model

The closed form problem is obtained by substituting equation (39) into equations (34) and (35) and further simplified. This leads to the volume averaged transport equation
(45)cαΦαεα∂tPαα+∇·εαcαqααPαα−∇·εαKα·∇Pαα+Kαη·∇Pβη−∇·εαuαdPβη−Pαα=∇·εαDα·∇Pαα+∇·εαDαη·∇Pηη+∇·uαcPηη−Pαα−SαηKαPηη−Pαα−εαuαs∇Pαα−εαuαβs∇Pηη−εαAvαηPηη−Pαα,where
(46)Kα=KαI+εα−1V∫AαηnαβbαdA,Kαη=εα−1KαV∫AαηnαηbαηdA,uαd=εα−1KαV∫AαηnαηsαdA,Dα=εα−1cαqα~bα,Dαη=εα−1cαqα~bαη,uαc=εα−1cαqα~sα,uαs=εα−1cαV∫Aαηnαη·bαqαdA,uαηs=εα−1cαV∫Aαηnαη·bαηqαdA,Avαη=εα−1cαV∫Aαηnαη·qαsαdA.Equation (45) is further simplified to obtain the macroscopic volume averaged transport equation for fluid flow in the α-region as follows:
(47)cαΦα∂tPαα+∇·qαePαα+uαηsPβη−∇·Kαe·∇Pαα+Kαηe·∇Pβη=∇·uαePβη−Pαα−AvαηPηη−Pαα−SαηKαPηη−Pαα,where
(48)Kαe=Kα+Dα,Kαβe=Kαβ+Dαβ,qαe=cαqαα+uαs,uαe=uαd+uαc,and the right-hand terms ∇·uαePηη−Pαα,AvαηPηη−Pαα and Sαηκα−1Pηη−Pαα represent the mass flux due to advective transport induced by α-η interaction at the interface, interfacial mass transfer by advection, and interfacial diffusive flux contribution at the interacting interfaces, respectively.Again, substitutingPη~ in equation (39) into equation (35) and simplifying leads to macroscopic volume-averaged transport equation in the η-region as follows:
(49)cηΦη∂tPηη+∇·qηePηη+uηαsPαα−∇·Kηe·∇Pηη+Kηαe·∇Pαα=∇·uηePαα−Pηη−AvηαPαα−Pηη−SηαKηPαα−Pηη,with similar representations of the coefficients defined for the α-phase equation. Due to the flow interaction between the α-phase and both the β- and γ-phases, it is important to account for flow interaction effects including cross diffusion and advection effects induced by both β-phase and γ-phase on the α-phase. The cross terms identified in the new equation (47) include uαηsPηη,Kαηe·∇Pηη, and ∇·uαePηη−Pαα where η∈β,γ and other terms arising from flow interaction between the α and both the β and γ phases. Thus, we obtain a system of volume-averaged transport equation for the α,β,γ-systems given as follows:
(50)cαΦα∂tPαα+∇·qαePαα+Uβγ−∇·Kαe·∇Pαα+Dβγ=∇·uαePββ−Pαα+∇·uαePγγ−Pαα−RαβPββ−Pαα−RαγPγγ−Pαα,(51)cβΦβ∂tPββ+∇·qβePββ+uβαsPαα−∇·Kβe·∇Pββ+Kβαe·∇Pαα=∇·uβePαα−Pββ−RβαPαα−Pββ,(52)cγΦγ∂tPγγ+∇·qγePγγ+uγαsPαα−∇·Kγe·∇Pγγ+Kγαe·∇Pαα=∇·uγePαα−Pγγ−RγαPαα−Pγγ,where
(53)Uβγ=uαβsPββ+uαγsPγγ,Dβγ=Kαβe·∇Pββ+Kαγe·∇Pγγ,Rαβ=−Avαβ−SαβKα,Rβα=−Avβα−SβαKβ,Rαγ=−Avαγ−SαγKα,Rγα=−Avγα−SγαKγ.The parameters obtained in this model and the associated parameter values used for the numerical simulation of the model are outlined in Table1. To accurately represent the permeability of the domain, the fracture and matrix permeabilities were evaluated from the relation proposed by [22], which are defined, respectively, as
(54)Kf=0.05758ϕf×hfeVf3f∈β,γ,K=ϕ42s2,where s is the pore area.Table 1
Table of parameters used in the model and parameter values for simulation [12, 23].
ParameterSymbolValue(s)Time (yrs)T5Fluid densityρ1Fluid viscosityμ0.001Porosity of domainϕϕα=0.40, ϕβ=0.05, ϕ=0.04CompressibilityCCα=0.15, Cβ=0.18, Cγ=0.20Fracture apertureeVeVβ=20, eVγ=55Fracture heighthfhβ=100, hγ=115


## 6. Results and Discussion

The VATE model obtained has new properties that makes it different from the existing advection/diffusion transport equations usually used for fractured porous media modeling. The cross diffusive terms with coefficientsKαβe,Kβαe, and Kαγe,Kγαe defined in this paper as
(55)Kije=εi−1KiV∫AijnijbijdA+εi−1ciqi~biji,j∈α,β,γ,i≠j,takes into account cross diffusive effect induced by one phase on the other. It is a physical flow property that is common with heterogeneous media which contributes to amount of substance transferred from one domain to another. In addition, cross-advection transport is another flow property that affects motion of fluid in heterogeneous media. In this new model, a cross-advective term with coefficient
(56)uijs=εi−1ciV∫Aijnij·bijqijdAi,j∈α,β,γ,i≠j,is obtained accounting for advective transport induced by one domain on the other interacting domain. One other important parameter obtained in this study is the diffusive coefficient which relates the mass of a substance that diffuses through a surface in a unit time to the concentration gradient of the substance. The effective diffusive coefficient
(57)Kie=KiI+εi−1V∫AijnijbidA+εi−1ciqi~bii∈α,β,γ,i≠j,obtained in this new model is significantly different from the diffusive coefficient in the existing models for multiscale modeling of fractured reservoirs. In this model, Kie,i∈α,β,γ captures both physical and geometrical properties such as viscosity, permeability, fluid density, and fluid velocity. The diffusive coefficient also takes into account advective and diffusive transport properties and may have either positive or negative values depending on the direction of the unit normal vector which determines the direction of flow of fluid.One challenging issue identified with multiscale fractured porous media modeling by multicontinuum approach is the design of the transfer function [4]. The existing and most used mass transfer function for multicontinuum modeling is that proposed by Barenblatt et al. [3] and Warren and Root [2], which according to Berre et al. [4], is less accurate for asymmetric and complex geometries or medium. In this new model, the resulting mass transfer term with coefficient
(58)Rij=−εi−1ciV∫Aijnij·qisidA−Sijεi−1Ki,takes into account both physical and geometric properties such as fluid viscosity, matrix and fracture permeabilities, flow velocity, fluid density, and porosity. The transfer term takes into account mass flow rate through the interacting interface from domain properties. Advective transport at the interface is as well captured in the transfer term.


### 6.1. Pressure Profile

The macroscopic volume averaged transport model for theα,β,γ-systems given in equations (50)–(52) was solved using the finite element method. To do this, the strong form of the problem (equations (50)–(52)) was first put in the weak form via variational formulation where a test function was chosen to multiply the strong problem. After descretizing the domain of interest, the unknown variable (pressure profile) was approximated using Galerkin approximation which led to a 3N×3N-system of equations and the solution of the problem sought for in the Hilbert-Sobolev space. The model was then simulated over a 5-year period in 0.6 monthly (about 2 weeks) interval, and the results obtained are shown in Figures 4–6. In Figure 4, the matrix pressure profile is presented showing a wave propagation from the inlet towards the outlet which turns radial behavior with higher pressure concentration at the center.Figure 4
Pressure distribution in the matrix domain with increasing time. A quasi-steady state regime observed after the first 2.4-month period.
(a)
Solution at 0.6 months(b)
Solution at 1.2 months(c)
Solution at 1.8 months(d)
Solution at 2.4 months(e)
Solution at 3.0 months(f)
Solution at 3.4 monthsFigure 5
Pressure distribution in the microfracture domain with increasing time. A quasi-steady state regime observed after 2.4-month period.
(a)
Solution at 0.6 months(b)
Solution at 1.2 months(c)
Solution at 1.8 months(d)
Solution at 2.4 months(e)
Solution at 3.0 months(f)
Solution at 3.4 monthsFigure 6
Pressure distribution in the macrofracture domain with increasing time. A quasi-steady state regime observed after 2.4-month period.
(a)
Solution at 0.6 months(b)
Solution at 1.2 months(c)
Solution at 1.8 months(d)
Solution at 2.4 months(e)
Solution at 3.0 months(f)
Solution at 3.4 monthsThis allows flow of fluid in two-dimensional,x,y-directions, allowing flow from the matrix into the surrounding fractures and also towards the outlet. A shorter transient flow behavior is observed after the first 2.4-month period followed by a quasi-steady state flow regime in time (as shown in Figure 7) which is typical of fluid transport in fractured porous media [24].Figure 7
Matrix pressure profile in increasing time showing an early quasi-steady state flow regime.The distribution of pressure in the fractures is presented in Figures5 and 6, both of which show horizontal-linear flow behavior following the direction of flow from the inlet (x=0) to the outlet (x=1) of the domain observed in Figure 8.Figure 8
Pressure profile showing the direction of flow of matrix, microfracture, and macrofracture domains.During the transient regime, all production is seen emanating from the fractures leading to a drop in pressure values in fractures. Thus, comparing matrix and fractures, fractures will enhance fluid production in fractured porous medium. The higher matrix pressure concentration value than fractures is as result of low permeable nature of matrix compared to fractures and their relative pore sizes.


## 6.1. Pressure Profile

The macroscopic volume averaged transport model for theα,β,γ-systems given in equations (50)–(52) was solved using the finite element method. To do this, the strong form of the problem (equations (50)–(52)) was first put in the weak form via variational formulation where a test function was chosen to multiply the strong problem. After descretizing the domain of interest, the unknown variable (pressure profile) was approximated using Galerkin approximation which led to a 3N×3N-system of equations and the solution of the problem sought for in the Hilbert-Sobolev space. The model was then simulated over a 5-year period in 0.6 monthly (about 2 weeks) interval, and the results obtained are shown in Figures 4–6. In Figure 4, the matrix pressure profile is presented showing a wave propagation from the inlet towards the outlet which turns radial behavior with higher pressure concentration at the center.Figure 4
Pressure distribution in the matrix domain with increasing time. A quasi-steady state regime observed after the first 2.4-month period.
(a)
Solution at 0.6 months(b)
Solution at 1.2 months(c)
Solution at 1.8 months(d)
Solution at 2.4 months(e)
Solution at 3.0 months(f)
Solution at 3.4 monthsFigure 5
Pressure distribution in the microfracture domain with increasing time. A quasi-steady state regime observed after 2.4-month period.
(a)
Solution at 0.6 months(b)
Solution at 1.2 months(c)
Solution at 1.8 months(d)
Solution at 2.4 months(e)
Solution at 3.0 months(f)
Solution at 3.4 monthsFigure 6
Pressure distribution in the macrofracture domain with increasing time. A quasi-steady state regime observed after 2.4-month period.
(a)
Solution at 0.6 months(b)
Solution at 1.2 months(c)
Solution at 1.8 months(d)
Solution at 2.4 months(e)
Solution at 3.0 months(f)
Solution at 3.4 monthsThis allows flow of fluid in two-dimensional,x,y-directions, allowing flow from the matrix into the surrounding fractures and also towards the outlet. A shorter transient flow behavior is observed after the first 2.4-month period followed by a quasi-steady state flow regime in time (as shown in Figure 7) which is typical of fluid transport in fractured porous media [24].Figure 7
Matrix pressure profile in increasing time showing an early quasi-steady state flow regime.The distribution of pressure in the fractures is presented in Figures5 and 6, both of which show horizontal-linear flow behavior following the direction of flow from the inlet (x=0) to the outlet (x=1) of the domain observed in Figure 8.Figure 8
Pressure profile showing the direction of flow of matrix, microfracture, and macrofracture domains.During the transient regime, all production is seen emanating from the fractures leading to a drop in pressure values in fractures. Thus, comparing matrix and fractures, fractures will enhance fluid production in fractured porous medium. The higher matrix pressure concentration value than fractures is as result of low permeable nature of matrix compared to fractures and their relative pore sizes.


## 7. Summary and Conclusions

This study has successfully modeled fluid transport by advection and diffusion in a three-phase multiscale fractured porous media consisting of a matrix and a 2-scale fractured network system with characteristic length scalesℓm and ℓM. A point advective-dispersive transport equation is obtained from the fundamental mass conservation equation together with the non-Darcy-Forchheimer momentum equation to incorporate into the model nonlinear pressure drop due to the high flow rate and permeability of the medium. The point advective-dispersive equation was spatially smoothed using a volume averaging technique, to obtain a macroscopic volume-averaged transport equation, from which a closure problem (equation (28)) was developed. A solution to the closure problem was proposed and written in terms of the averaged quantities and some closure variables that are chosen based on the boundary-value Problems 1–3. The resulting macroscopic equation has as follows:
(i)
Cross-advective and diffusive terms that account for induced fluxes between the interacting domains(ii)
Effective diffusive coefficients, which incorporates domain properties such as viscosity, fluid density, and flow velocity(iii)
A new representation of the mass transfer function, which has additional geometric and physical properties including flow velocity, compressibility, fluid density, and porosity. The area integral and the surface normal vectors deal with variable fracture surface area and fracture orientation, respectivelyThe results of the numerical simulation of the model exhibit the radial, horizontal-linear flow behavior and a transient and quasi-steady-state flow regime that is typical of naturally fractured porous media.



---

*Source: 5731988-2022-03-31.xml*