A new macro-scale volume averaged transport model for diffusive dominated non-Darcian flow problem in multi-scaled naturally fractured reservoirs

Diffusive transport in porous media is a complex process in multi-scaled fractured media modeling. This paper presents a diffusive transport model for non-Dacian flow in a naturally fractured reservoir with triple porosity and permeability. To address the non-Darcian flow behavior associated with fluid transport in fractured porous media, the Darcy/Forcheimer equation was used. A point-diffusive equation was obtained from mass conservation and the Darcy–Forcheimer momentum equation; this is used together with interface conditions to incorporate the microscopic properties of the domain. Subsequently, the resulting equation was spatially smoothed to obtain an effective macroscopic average model. The macroscopic model obtained, unlike the existing models, has a cross-diffusive term for mass transport by induced fluxes and a mass transfer term accounting for mass transfer between the matrix and the surrounding fractures via the interface. The numerical simulation displayed a horizontal-linear flow behavior in the fractured network instead of a radial flow in the matrix. The results further suggest that despite the fractures aiding in fluid transport, they enhance fluid production in the reservoir compared to the matrix.


Introduction
In recent times, modeling fluid transport in fractured reservoirs piqued the interest of oilfield experts and professionals in reservoir modeling. Fractures existing in reservoirs have a high-level effect on transport processes in reservoirs. It introduces secondary conduits, which makes the porous media more permeable, and hence more productive (Zhao et al. 2022). The fractures, on the other hand, have an adverse effect of increasing the rate of decline of fluid production, which may lead to early breakthrougha phenomenon in petroleum production that requires early injection (Bratton et al. 2006). The complexity of modeling fluid transport in fractured reservoirs increases with the increasing level of fractures and their interaction with the matrix. The geometrical and physical properties of the fractured domain influenced by the presence of fractures also poses a significant challenge in reservoir modeling for diffusive transport processes. The inception of fractured reservoir modeling is identified in the works of Barenblatt et al. (1960) and Warren and Root (1963), which considered modeling transport processes in the fractured medium as a dual continuum for the rock/matrix and fractures, with each continuum exhibiting different flow properties. This approach led to the development of dual-porosity and dualpermeability models. Consequently, the triple-porosity and triple-permeability models for multi-scale fractured Isaac Kwame Dontwi and Peter Amoako-Yirenkyi contributed equally to this work.
1 3 reservoirs with the fractures considered at different lengthscales, and depending on fracture size were also developed.
The existing fractured reservoir models are based on laminar flow regime (Choquet 2004;Nakshatrala et al. 2016;Wei et al. 2018;Zhang et al. 2018), with few considering dual-porosity models for nonlinear flows (Altinörs and Onder 2010). According to Altinörs and Onder (2010) and Berre et al. (2019), the use of laminar flow assumption for fractured reservoir modeling is inaccurate, due to the high flow velocity and permeability in fractures, resulting in nonlinear flow behavior in fractured reservoirs (Barros-Galvis et al. 2018). Another limitation identified with existing models is the representation of the transfer term, which accounts for mass transfer between the interacting domains; the matrix and its surrounding fractures. In existing models, mass transfer terms like those proposed in the work of Barenblatt et al. (1960) and Warren and Root (1963) for multi-continuum modeling only accounts for geometrical properties such as fracture-matrix interface area and physical properties like fracture-matrix permeability and fluid viscosity (Berre et al. 2019). Moreover, in reactive systems, the gradient in concentrations in each of the interacting media induces a flux in the domain at the matrix-fracture interface. As a result, mass transfer between interacting domains may be affected not only by the geometrical and physical properties of the fractured domain, as well as the distribution and orientation of fractures, as reported in Berre et al., (2019) but also by cross-diffusion effects due to induced fluxes. As a result, in order to accurately represent the mass transfer term, the transfer coefficient must be extended or generalized to capture additional important domain properties. Diffusive transport in multi-phase systems has been extensively studied using the volume averaging theory concept (Borges da Silva et al. 2007;Whitaker 1999), which has the property of up-scaling by taking micro-scale information from the multi-phase system to the macroscopic scale level (Francisco et al. 2017). The method of volume averaging was successfully applied to transport problems in heterogeneous porous media for active and passive dispersion n the works of Ahmadi et al. (1998), Carbonell and Whitaker (1983) and Whitaker (1999). This technique led to volume-averaged transport equations with effective diffusive coefficients and mass transfer terms which take into account additional geometrical and physical properties of the domain.
In this work, fluid transport process for nonlinear flows in a multi-scale naturally fractured reservoir where the medium is assumed to have two different length-scales of fractures, m (micro-fracture) and M (macro-fracture), is considered. The matrix and the two scales of fracture networks are assumed to have different flow properties, in this case, thus making the medium heterogeneous. Unlike the existing models, the proposed model began with the development of 'point' advective diffusive transport equation from the mass conservation equation and the Darcy-Forchheimer momentum equation in order to address the issue of nonlinear flow behavior in fractured domains. The volume averaging technique was then applied to the point-diffusive transport equation at the microscopic level to obtain a novel local volume averaged triple porosity and permeability model representing diffusive transport at the macroscopic level. The resulting model incorporates non-Darcian behavior by using the Darcy-Forchheimer momentum equations. It also has cross-diffusive terms that take into account crossdiffusion between interacting phases of the medium and a mass transfer term that considers the domain's diffusive properties. Flow properties such as flow velocity and fluid viscosity are well captured in the transfer term, and it also addresses the issue of fracture orientation and varying fracture interface areas. A detailed formulation, analysis and numerical results of the developed model are presented in the subsequent sections.

Micro-scale boundary value problem for non-Darcian flow in multi-scale fractured porous media
We consider a multi-scale fractured porous media with three phase system, namely, the matrix, the micro-and the macrofracture networks; denoted as -, -and -phase, respectively, illustrated in Fig. 1.
The diffusive transport process, which involves the mass transfer between the matrix and the 2-scale fracture networks were investigated in this study. At the matrix-fracture interface, we proposed a sharp interface of zero thickness with no adsorption or desorption. The instantaneous flow velocities at the interacting interfaces ( − and − interfaces) equal the velocity normal to the interface. As a result, the pressure and the normal component of the pressure fluxes at the interacting interfaces are both continuous. The single-phased fluid transport processes with non-Darcian flow behavior in the three-phase system are modeled as where Φ, , and are, respectively, the porosity of the domain, fluid density and flow velocity. Fluid transfer occurs in fractures with a flow rate Q, that is proportional to the pressure drop, mathematically represented as; Darcy empirically established an equation relating the flow velocity and the pressure drop given in Eq. 3 (Darcy et al. 1856) (1) where K and are, respectively, domain permeability and fluid viscosity. The concept of permeability in porous media has been empirically proven to be valid in two important conditions, thus; 1) the fluid must be Newtonian, and 2) inertial forces must be small enough in comparison to viscous forces. This proves that for a porous medium with a characteristic pore scale l, the Reynolds number, must be small enough where , v, and represent the fluid density, characteristic velocity and viscosity, respectively. As a result of the viscous drag forces which may be dominant in the porous medium, the condition for the validity of Darcy flow in a porous medium is According to Dukhan and Minjeur (2011), porous media exhibits different permeabilities depending on the flow velocities of the fluid in the medium. Hence, we may have low flow velocities accounting for low permeable medium while high flow velocities gives rise to high permeable porous medium. Therefore, when a fluid in motion transits from low permeable medium to high permeable medium, the pressure drop ceases from being linear to nonlinear. In this flow regime, the Darcy flow concept no longer holds for high Re (Zolotukhin and Gayubov 2022). As a result, we use the Darcy-Forcheimer equation (5); which incorporates the inertial term, ‖ ‖ that accounts for the nonlinear pressure drop where = c K a Φ b is a function of porosity and permeability, K a is inertial permeability, and is the flow velocity. Using the mass conservation Eq. 1 and the Darcy-Forchheimer equation 5, we obtain the diffusive transport equation for non-Darcian flow processes as; where the diffusivity coefficient resulting from the nonlinear flow behavior is obtained as The governing equations associated with the diffusive transport process considered in the three-phase system described in Fig In addition to the interfacial boundary conditions BC1 and BC2, BC3 and BC4 are the boundary conditions at the exits and entrances of the macroscopic region. It is important to note that BC3 and BC4 are only known in terms of average quantities at the macroscopic level. A e , A e denote the area of the entrances and exits of -region and -region, respectively, at the boundary of the macroscopic domain.

Spatial smoothing
To obtain the macroscopic transport equations that capture microscopic properties of the domain at the macroscopic level, the volume averaging method is employed. The macroscale variables are defined through spatial smoothing (averaging) of the domain from the microscopic level. Here, we begin the averaging process by associating with every point in the macroscopic region an averaging volume (illustrated in Fig. 1) that is constant in space and time. Given the threephase system, we express the averaging volume, V.
Thus, the volume fraction, i for each phase i in the averaging volume can be expressed as and further express the superficial average and intrinsic average of the pressure distribution in the averaging volume, respectively, as; these are related by Eq. 17.
The average of the spatial derivative is defined using the spatial averaging theorem as To begin with the up-scaling process, we consider first the flow interactions between the -region and the -region. Following the volume averaging theory and since the averaging volume, V , is invariant in time, the average of the mass accumulation term can be simplified as The spatial averaging theorem in Eq. 18 can be used to obtain the average of the diffusive term as; which subsequently simplifies to Using Eqs. 19 and 21 with similar simplifications for the -phase, the spatially smoothed equations for flow interaction between the -region and the -region are given in the form (in the −region) and in the -region we have where the interfacial flux, ing to mass transfer from -region into -region, is represented as (Quintard and Whitaker 1993) to connect the ( , ) regions at the interface. The S v in the coefficient is the total interfacial surface area at the ( , ) -interface. In this representation (Eqs. 22 and 23), the preferred intrinsic pressure variables describing the diffusive transport process is used, and the spatial deviation variables P and P are used in the diffusive integral terms according to the relation The spatial deviation variables present in the volume-averaged Eqs. (22) and (23) are to be estimated in terms of the averaged quantities, by developing a closure problem. This is done to obtain a closed form averaged transport equations at the macroscopic region.

Closure problem
In order to derive a macroscopic transport equation written only in terms of the averaged variables, we develop a closed form problem of the average transport Eqs. (22) and (23). This is made possible by developing a closure problem from which P = ⟨P ⟩ +P and P = ⟨P ⟩ +P .
the deviation variables can be estimated in terms of the averaged variables. By dividing the average transport equation by the phase volume fractions and ; and subtracting it from the original (point) equation, the resulting partial differential equation for the spatial deviation variable is derived in the phase as; with a similar representation in the -phase.
To solve the spatial deviation equation at the microscopic level, we identify and work with terms that contribute to the diffusive transport equation, thus 'non-contributing' terms are removed from the equation by estimating and comparing terms at length-scales (microscopic scale) and (macroscopic scale), of the measuring volumes. From this, we obtain the boundary value problem for the spatial deviation variables as (in the -phase); and in the -phase we have Obviously, the solution to the boundary value problem for the deviation variables depends on the source and the boundary data which are represented in the problem in terms of ∇⟨P ⟩ , ∇⟨P ⟩ and (⟨P ⟩ − ⟨P ⟩ ).
Following the works of Whitaker (1999), Carbonell and Whitaker (1984) and Ahmadi et al. (1998), the spatial deviations are expressed as functions of the macroscopic variables in terms of the closure variables and s, written as;

The closure variables
The closure variables s i and i , ij (i ≠ j), i, j ∈ { , } are, respectively, scalar and vector fields, that, respectively, maps the averaged quantities ⟨P i ⟩ i and its spatial derivative ∇⟨P i ⟩ i unto P i which is specified according to the boundary value problems (I-III), defined for the closure variables from the closure problems 27 and 28. The boundary value problem is developed by substituting the proposed solution given by Eq. 29 into the closure problem represented by Eqs. 27 and 28. From this, we obtain the boundary value problem I as;

Problem I
Here, the boundary conditions BC3 and BC4 in Eq. (28) are replaced by spatially periodic conditions which allows us to write the periodic condition (Periodicity) for the closure variables in Eq. 31. The periodicity condition is imposed since the boundary problem for the spatial deviation variables is to be solved for P , at the microscopic scale of characteristic length , and not the entire macroscopic region. To do this, we abandoned BC3 and BC4 in Eq. (28) imposed at A e since it influences the P -field in the representative microscopic region on the order of and replaced it by spatially periodic condition, if the representation is treated as a unit cell in spatially periodic porous medium (Whitaker 1999). Similar simplification is applied to obtain the second and third boundary value problems for the closure variables as; Problem II (29)

The closed form problem
The closed form problem representing diffusive transport equations at the macroscopic region is obtained by substituting the proposed solution in Eq. (29) into the volumeaveraged Eqs. (22) and (23). The diffusive term with the deviation variable can be simplified as follows; where Substituting the diffusive term into Eq. (29) and the coefficient of the mass transfer term, R is obtained as; Similarly, the -phase macroscopic transport equation is obtained as; Next, spatial smoothing for transport equations for flow interactions between the -phase and the -phase also leads to -phase and -phase macroscopic transport equations, respectively, obtained as (in the -phase); and in the -phase as; where the boundary value problem for the closure variables s i , i , and ij (i ≠ j) i, j ∈ { , } are defined similarly to problems I-III obtained in "The closure variables" section.

The macro-scale transport equations for multi-scale fractured reservoirs
To fully represent transport processes in the model proposed for the , , -system, the two cross-diffusive terms and the mass transfer terms (i.e., ⋅ ∇⟨P ⟩ , ⋅ ∇⟨P ⟩ and R (⟨P ⟩ − ⟨P ⟩ ) , R (⟨P ⟩ − ⟨P ⟩ ) ) appearing in the volume averaged transport equation in the -region for flow interaction between the and that of -and -region needs to be accounted for in the matrix model. This leads to a system of volume averaged transport equation for diffusive transport in the multi-scale fractured reservoir as; where the effective diffusivity and cross-diffusive coefficients are, respectively, represented as; (for i, j ∈ { , , } , i ≠ j) and the coefficient of the mass transfer term, R ij is represented as;

Input parameters
The input parameters for the numerical simulation are shown in Table 1. To accurately represent the permeability of the domain for the numerical simulation, these values were evaluated from the properties of each of the domain using (45) the relations proposed in the work of Tiab and Donaldson (2016) for fractures and matrix, respectively, defined as where the micro-and macro-fracture heights are given, respectively, as h = 100 and h = 115 , and the fracture aperture e Vf defined in Table 1.

The cross-diffusive terms and transfer function
The cross-diffusive term in the averaged transport equations results from flux contribution induce by one domain on the other interacting domain. Since the concentration gradients in either of the interacting domains induces a flux on each other, the cross-diffusive term becomes an important component in this model, which describes diffusion between the matrix and the fractures. It is observed that the crossdiffusive term with coefficient contributes to the fluxes of each interacting domain (describing the diffusion between matrix and fracture) by a significant amount when compared to the effective diffusive term with coefficient The effective diffusive coefficient obtained in this model has surface normal component ij , that depends on flow velocity, fluid density and viscosity and porosity of domain, which significantly affects pressure distribution of the domain.
The transfer term with coefficient obtained in this model, captures, in addition to the geometrical and physical properties of the domain, diffusive properties at the interacting interface, which have high influence on fluid transport processes, and the surface normal vector that deals with issues of fracture orientation. Flow velocities, i of the interacting domains is another flow property captured in the transfer function which affects mass transfer process by either increasing or decreasing the diffusive coefficient D i . The transfer function obtained in this model is more accurate compared to that proposed by Barenblatt et al. (1960), which is less accurate according to Berre et al. (2019).

Pressure profiles
We discuss in this section the pressure profiles from the numerical simulation of the proposed model. A fractured reservoir producing at a constant rate is considered, where the injector and the producer are set at x = 0 and x = 1 , respectively. The numerical results were obtained using the finite element method via variational formulation of the strong form of the problem and the solution sought for in a Hilbert-Sobolev space. Figure 2 shows the flow distribution in the matrix domain of the fractured porous medium. The matrix distribution shows radial movement from a higher concentrated region toward the surrounding fractures and/ or the exits of the domain. It is observed that, with time, the pressure distribution in the matrix domain exhibits a slightly short transient behavior preceding a quasi-steadystate regime-a phenomena which is identical to that of transient fractured porous media models (Bourdet et al. 1984;Warren and Root 1963). The flow transport process in the fractures for micro-and macro-fracture domains is shown, respectively, in Fig. 3 and 4, exhibiting radial linear flow movement for both domains. Fractures act as conduits and aid in fluid production. A higher concentration of the pressure values is observed from the pressure profiles at the inlet (the injector, x = 0 ) with decreasing concentration across the flow path to the outlet (the producer, x = 1 ). This allows mass transfer from the inlet to the outlet of the domain in a radial, linear behavior. Due to the radial, linear movement in the fractures compared to the radial movement in the matrix, mass transport is also high, enhancing productivity. We note that a steady-state flow regime is also observed in the fractures after a short transient flow regime.

Conclusions
Fractures in porous media occur at different scales, and this introduces high heterogeneity due to the different flow properties that occur at different regions of the fractured domain. Moreover, due to the presence of fractures, high flow velocity and permeability occur, causing a deviation from the laminar flow regime to non-Darcian flow behavior. As a result, we propose a multi-scale fractured reservoir model to study the complex transport process in naturally fractured reservoirs, where non-Darcian flow behavior is captured using the Darcy-Forcheimer flow equation. The model is a multi-continuum model consisting of the matrix surrounded by micro-fracture and macro-fracture networks leading to triple porosity and triple permeability model. In this model, each continuum has its own porosity, permeability and compressibility properties. An up-scaling technique is used to accurately represent diffusive transport processes and mass transfer between interacting phases for the proposed fractured reservoir model, which captures flow properties from the microscopic level to the macroscopic level. The local volume-averaged transport equations (VATE) proposed in our study have better features and are more robust for studying diffusive transport processes in naturally fractured reservoirs than the existing triple-porosity permeability model like that proposed in Whitaker (1999), Wu et al.  (2004) and Zhang et al. (2018). These existing models has no cross-diffusive properties, has constant diffusive coefficients that are independent of microscopic properties, including fluid density, flow velocity and fluid viscosity, and also assume the flow to be linear, which according to Bere et al., (2019) is inaccurate. Numerical simulations of the proposed model capture the expected behavior of pressure profiles in multi-scaled naturally fractured reservoirs.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.