An integrated approach for relative permeability estimation of fractured porous media: laboratory and numerical simulation studies

Most carbonate reservoirs in Middle East are characterized as porous fractured reservoirs. Estimation of relative permeability of these highly heterogeneous reservoirs is challenging due to the existence of discontinuity in the fluid flow fractured porous media. Although relative permeability is an essential data for simulation of flow in fractured media, few attempts have so far been made to estimate the relative permeability curves. Most notable are the studies by Akin (J Pet Sci Eng 30(1)1–14, 2001), Al-sumaiti and Kazemi (2012), and Fahad (2013). This paper presents an integrated approach to history matching the oil drainage tests, which were carried out by unsteady state, on glass bed models with a single fracture at different orientations and to estimate the relative permeability curve. The integrated approach includes an inversion algorithm coupled with forward numerical modeling of fluid flow. The history matching of the displacement test data was obtained by using the Levenberg–Marquardt algorithm to minimize the error between the simulated and experimental data. In this algorithm, Corey-type power law is used to create relative permeability curves during the optimization procedures. The forward modeling is a 3D multiphase fluid simulator for flow through discrete fractures. Numerical results of fluid flow profiles and the optimized relative permeability curves for single fracture with different orientations and experimental validation with oil drainage tests are presented. The results of the optimized relative permeability data for single fracture are in a good agreement with the data derived by the correlation of Fahad (2013). These results prove that the presented approach can be used to upscale the relative permeability curve from laboratory scale to reservoir grid scale. The work on the upscaling of the estimated relative permeability curve of fractured porous media is under preparation and will be published soon.


Introduction
Relative permeability plays a critical role in evaluating the production potential and recovery factor for conventional and unconventional reservoirs, and yet, it is considered to be one of the most uncertain parameters in the governing equations of multiphase fluid flow. There are two methods usually used in determining relative permeability in the laboratory by fluid displacement experiments on core samples: (a) in steady-state and (b) unsteady-state conditions. The laboratory-derived relative permeability by steady or unsteady state is upscaled to reservoir scale by history matching of the blocked-based simulated pressure and production data with that recorded during production. In a fractured reservoir, the evaluation of relative permeability is challenging due to the existence of discontinuity in the fluid flow through fractured porous media. The linear relative permeability model (analytical) proposed by Romm (1966) has been used in fluid flow simulation in fractured reservoirs. In this approach, different forces, such as the capillary, viscous, and gravity, are ignored which lead to significant errors in the prediction of fluid recovery (De la Porte et al. 2005) particularly in fractured and unconventional reservoirs. Recent studies have shown that the relative permeability in fractured reservoirs is not a linear function of saturation (Chima et al. 2010;Fahad 2013).
The previous experimental work that has been carried out on fractured cores can be classified into two groups: (a) experiments based on the assumption that the fluid flow is dominant in fractures and (b) experiments based on the assumption that the fluid flow is dominant in porous matrix and fractures. Romm (1966) presented an analytical equation to calculate oil-water relative permeability in fractures based on experimental results using parallel glass plates with smooth surfaces. Based on displacement tests of using two parallel glass plates (kept in vertical position) Romm (1966) concluded that the relative permeability is a linear function of saturation. It was also noted that the linear relative permeability curve is not applicable for wellconnected fractured system. Pruess and Tsang (1990) presented a conceptual numerical model for multiphase fluid flow through rough-walled fractures. The model utilized the quantitative description of the fracture pore space in terms of aperture distribution, which was obtained by direct measurements on fracture surface. The void space of a rough-walled fractures is then conceptualized as a twodimensional heterogeneous porous media. This study showed that the relative permeability strongly depends on fracture aperture. In addition, the authors concluded that the sum of wetting-and nonwetting-phase relative permeabilities within the fracture is less than one at intermediate saturation due to the strong interference between fluid phases. Persoff et al. (1991) studied two-phase flow for gas and liquid in a rough-walled fracture with controlled injection rates and pressures using Hassler and Brunner (1945) method and porous end caps to avoid capillary end effects and to maintain uniform capillary pressure across the entire fracture plane. The results show that there is strong phase interference and the relative permeability is not a linear function of saturation. Pieters and Graves (1994) used a video imaging technique to measure saturations in rough-walled rock fracture, and then, the twophase relative permeability is calculated based on Welge's (1952) interpretation of the Buckley-Leverett theory after calculations of the fractional flow. Image processing was performed at various time intervals during the displacement process (before and after breakthrough time). The authors concluded that the relative permeability is not a linear function of saturation. Chima et al. (2010) derived an analytical equation to calculate the relative permeability curve in fractured porous media for two-phase fluid flow. The authors used the mass and momentum balance equation for isothermal system, Newton's law of viscosity, and Darcy's law in the derivation of the proposed equation. The proposed equation related the relative permeability to saturation and fluid viscosity. The authors observed that the analytical equation gave results close to that of Pieters and Graves (1994) and the phase relative permeability in fractured porous media is not a linear function of saturation.
The notable two-phase displacement experiments that are conducted in fractured porous media include the work of Akin (2001) andFahad (2013). Akin (2001) used a history matching technique to calculate relative permeability curve from steady-state experimental data which was carried out on a core sample with a single vertical fracture. The author used a dual-porosity/dual-permeability approach to simulate fluid flow. Akin (2001) found that the fracture relative permeability is a power law function of the corresponding phase saturation. Fahad (2013) performed unsteady-state drainage tests on a glass bead model using single and multiple fractures (both parallel and intersected) with different orientation angles. Based on the results, the authors proposed a new equation for relative permeability for fractured porous media. The equation relates the relative permeability to the number of fractures and fracture orientations. The experimental results also showed that the relative permeability is not a linear function of phase saturation.
Three-dimensional simulation of multiphase fluid flow in fractured porous media Simulation of multiphase fluid flow in naturally fractured reservoirs is a challenge due to complexity associated with flow processes in porous matrix fracture system. A detailed study of multiphase fluid flow for estimation of production potential of fractured reservoirs is not well documented. In this paper a full derivation of the multiphase fluid flow equations in a poroelastic environment is presented to accurately reflect the fluid flow behavior and evaluate fluid recovery under different driving mechanisms. A number of different approaches have been used to simulate immiscible multiphase fluid flow in fractured reservoirs. Among these approaches most notable is the dual-porosity/dualpermeability one.
In dual-porosity/dual-permeability approach, Kazemi et al. (1976) developed a three-dimensional numerical simulator for water/oil system in fractured reservoirs. The equations used to simulate two-phase flow are extensions of Warren and Root's (1963) single-phase flow. The numerical model accounted for gravity, viscous, and capillary forces. Governing equations for matrix porous media and fractures are solved separately by Gauss-Seidel iterative technique using finite difference method. In this approach, the fracture pressure and water saturations in fracture and the matrix pressure and water saturation in matrix are solved iteratively. This process is repeated until a certain convergence is achieved. Kazemi et al. (1976) did not included the gravity term, and the same capillary pressure curve was used for both fracture and matrix systems. Thomas et al. (1983) extended the approach that has been presented by Kazemi et al. (1976) by replacing the two-phase fluid flow with three-phase flow in the derived equations. In their approach the authors assumed that the flow between the matrix and fractures is governed by a local transfer function and that there are no communications between matrix blocks. The mathematical formulations were solved implicitly for pressure, water saturation, and gas saturation. The authors found that the implicit treatment of matrix/fracture flow resulted in a more stable and efficient model compared to the method proposed by Kazemi et al. (1976). In addition, Thomas et al. (1983) concluded that the gas saturation has a significant impact on oil recovery when used for water flooding. Lewis and Ghafouri (1997) extended Thomas et al. (1983) model by accounting for poroelasticity. The mathematical equations are derived based on geomechanics equilibrium and multiphase mass conversation. The authors discretized the equations using Galerkin finite element formulations. Lewis and Ghafouri (1997) found that the coupling effect has a significant impact on the reservoir fluid flow behavior.
In conclusion, all the above-mentioned models have used dual-porosity/dual-permeability approach in simulating multiphase fluid flow through naturally fractured reservoirs. These models, however, are not physically and mathematically rigorous for multiphase fluid flow as the derived equation is direct extensions of single-phase flow equations and ignored the mechanisms that govern the multiphase fluid flow. In other words, multiphase transfer function used in dual-porosity/dual-permeability approach is a direct generalization of single-phase transfer function and the effect of capillary force is neglected. Karimi-Fard et al. (2004) simulated the two-phase fluid flow in fractured reservoirs by integrating finite volume and Galerkin finite element methods to guarantee mass conservative and provide direct physical interpretation of the discretized governing equations. In order to obtain a continuous pressure and saturation profile, Hoteit and Firoozabadi (2005) developed a compositional simulator based on mixed finite element and discontinuous Galerkin method. Firoozabadi (2004, 2007) used the capillary continuity equation and capillary pressure saturation function at fracture matrix interface to calculate the saturation in fracture and matrix blocks within the control volume. Raviart and Thomas (1977) and Arnold et al. (2002) used mixed finite element method (MFE) to approximate velocity fields in highly heterogeneous porous media. This method has been widely used to simulate single-phase fluid flow through fractured porous media (Vohralík 2007). Hoteit and Firoozabadi (2008) extended (MFE) to two-phase fluid flow in reservoirs with few discretely oriented fractures including gravity and capillarity effects. Their model accounted for fluid flow through matrix, fracture, and at matrix fracture interface. Monteagudo et al. (2011) and Moinfar et al. (2014) simulated multiphase fluid flow coupled with geomechanics in few discretely oriented fractures. Monteagudo et al. (2011) used finite volume method to solve the flow equations and Galerkin finite element for the geomechanics. Moinfar et al. (2014), on the other hand, proposed a fully implicit compositional reservoir simulator using embedded discrete fracture model (EDFM) to simulate complex displacement process (e.g., miscible gas injection) in fractured reservoirs. The authors used the empirical model proposed by Bandis et al. (1983) which is called ''nonlinear Barton-Band joint model'' to relate the effect of effective normal stress to the normal closure of the fracture. The problem with Barton-Band joint model is its underestimation of hydraulic aperture value (Piggott and Elsworth 1991). In all of these studies the authors simulated fluid flow in few fractures.
In summary, most of the numerical models used to simulate multiphase flow in fractured porous media are uncoupled with geomechanics problem or coupled but using dual-porosity/dual-permeability approach. Also, these models have different limitations of incorporating full permeability tensors and capable of simulating only few numbers of discretely oriented fractures. In addition, application for multiphase upscaling coupled with geomechanics problems has not been addressed.
In this paper, a fully coupled multiphase fluid flow in discretely oriented fractured porous media in a poroelastic frame work will be presented. The simulation workflow for two-phase fluid flow is based on upstream flux weighted finite element discretization method in a poroelastic framework. Moreover, the viscous and the gravity effects are considered in mass and momentum balance equation. The capillary pressure effect is included during the discretization of the partial differential equations of multiphase fluid flow. In flow simulation the capillary pressure saturation function is used for water saturation distribution calculations at fractures and matrix interface. The developed numerical scheme is validated against an analytical solution for homogenous porous media only.
Next, an inversion algorithm will be presented and then integration between the developed multiphase flow simulator and the inversion technique will be used to estimate relative permeability curve for laboratory drainage tests data (production data) that carried out on glass bed models (Fahad 2013). J Petrol Explor Prod Technol (2020) 10:2499-25162501 Governing equations for multiphase fluid flow in a poroelastic framework In general, behavior of two-phase fluid flow system through fractures and matrix porous media is governed by generalized Darcy's law and continuity equation (see ''Appendix'' for a comprehensive derivation). The Darcy's law is expressed as: General continuity equation for wetting phase incorporating the concept of effective stress can be expressed as follows: General continuity equation for nonwetting phase can be expressed as follows: where / is the porosity of the media, c is the rock compressibility, p is the average pressure, s p is the saturation for each phase, u ps is the relative velocity vector between fluid phase and solid phase, k ij is the permeability tensor, k rp is the relative permeability for each fluid phase p, l p ; q p , and p p are dynamic viscosity, density of fluid, and fluid pressure for each phase, respectively, g i is the gravity acceleration vector, b p is the fluid formation volume factor, K m is the bulk modulus of solid grain, D is the elastic stiffness matrix, and Q p represents external sources or sinks.

Saturation equations
Oil saturation equation: Water saturation equation: where s i and v i are the saturation and velocity for phase (i), respectively.
To obtain the numerical solution of this highly nonlinear system of equations, suitable initialization and boundary conditions are designated first, and then some of auxiliary functions are employed which are known as the constitutive relationships.

Auxiliary equations
The porous medium voids are assumed to be filled with water and oil, and the sum of their saturation is unity: The pressure in each phase in porous media is related though the capillary pressure equation as follows (assume that the media are water wet): where p cow is the capillary pressure between oil phase and water phase.

Validation of 3D poroelastic numerical model
Verification of the developed 3D two-phase numerical model is necessary to ensure that the three-dimensional concepts are implemented correctly. In view of this first the elastic and two-phase fluid flow problems are validated separately. For the validation of elastic problem, a complex test using 3D beam element is used, while the fluid flow problem is validated against an analytical solution of Buckley-Leverett.

Validation of elasticity problems using FEM tool box
A 3D beam element is used with dimensions of (0.1 m 9 1 m 9 0.1 m) to test whether the developed 3D elastic numerical model is capable of predicting the stresses profile for complex loading schemes. Shear traction of 450 N/m 2 is applied on one side along y direction, and a fixed displacement is applied across x-z plane (see Fig. 1). The Young's modulus and Poisson's ratio are 200 GPa and 0.2, respectively. Results of this study are compared with that from a FEM tool box and presented in Fig. 2a, b. From the figure, it can be seen that the displacement and stresses profile resulted from the elastic numerical model and FEM tool box are very close to error ranging from 5 to 10 %. This level of error can be attributed to the use of gauss points: Gauss points used in the numerical model and FEM tool box are 15 and 4, respectively, and also the different methods used for stresses calculation. Noteworthy, the elastic numerical model produced the highest stresses (around 135 Pa) along the top edge of the 3D beam, particularly near the corners which is 5 Pa higher than that from FEM tool box, while the lowest stresses (around -450 Pa) were produced on the bottom edge of the 3D beam which is -45 lower than that from the FEM tool box. In addition, the developed numerical elastic model gave a continuous stress profile better than that produced by FEM tool box.

Two-phase numerical model validation against 1D Buckley-Leverett solution
The two-phase fluid flow numerical model is validated against 1D Buckley-Leverett analytical solution (Buckley and Leverett 1942) using simple 1D homogenous reservoir. The reservoir domain is divided into a number of elements with one producer and one injector. The injection rate is set to 1000 bbl/day and the initial reservoir pressure to 5000 psi, while the producer is producing at a variable rate at a constant bottomhole flowing pressure of 500 psi. The viscosities of wetting phase and nonwetting phase are 1 cp and 0.7 cp, respectively. In this case the capillary pressure effect is ignored; therefore, an explicit equation for water saturation calculation has to be used. Discretization of water saturation equation using standard finite element method produces a solution with spatial oscillations due to its hyperbolic nature. To overcome this, Galerkin least square technique (GLS) is employed to stabilize the solutions. Figure 3 shows the water saturation profile for Buckley-Leverett and finite element numerical solutions using different numbers of mesh elements. As can be seen from Fig. 3, with the increasing of mesh elements number, the water saturation profile becomes close to the analytical solution. In addition, the numerical saturation flood front does not exhibit the same piston-like displacement that is shown by the analytical solution. This is mainly because of the numerical dispersion due to time and space discretization. Numerical dispersion depends on element size, time step size, velocity of frontal advance, and numerical formulation. In addition, in finite element method, shape functions are used to approximate the field variables (pressure and saturation). Since these shape functions cannot capture the field variation exactly, the resulting interpolation error leads to numerical dispersion. Figure 4 shows the stable water saturation profile at different points away from the injector with time for 300 elements mesh.

Relative permeability measurement using a glass bead model
In order to improve our understanding of the two-phase flow process in oriented and intersected fractures, an innovative laboratory methodology was designed and developed by Fahad (2013) using a glass bead models of size 20 cm 9 10 cm 9 0.2 cm. The fractures in these models are made using two parallel glass plates separated by microfine glass beads of 2 mm width. In these packed microfine beads, single and multiple fractures of 0.1 mm 9 2 mm 9 100 mm are made by large diameter glass beads. Fluids are injected at a constant injection rates by using a piston pump. A digital camera was used in order to take pictures of glass bead model during the displacement process. A schematic representation of the experimental setup is presented in Fig. 5. Water and Soltrol-130 with 1.4 cp and 1.002 cp, respectively, are used during the unsteady-state displacement test. Constant oil injection rates of 1 and 2 cm 3 /min were used to displace the water phase, and the volume of produced fluids was recorded with time. These procedures were repeated by using different fracture system, such as multiple fractures with variable orientation and intersections. With the use of analytical equations introduced by Johnson, Bossler and Neumann (JBN) method (1959), the authors were able to calculate the relative permeability for each fracture system. Finally, the authors came up with a new correlation to relate the relative permeability of oil and water with number of fractures and their orientations and intersection.
The drainage test is carried out on a single and multiple fracture systems as shown in Fig. 6. The authors carried out a series of tests using single and multiple fractures with different orientation angles (see Fig. 6), and further details can found in the original literature. In this study a single vertical fracture with 0°in the direction of flow (see  data for different fracture systems. For this purpose, a history matching optimization algorithm that incorporates the multiphase numerical model presented in this paper. Moreover, in this section Fahad (2013) relative permeability correlation is evaluated with a view that the correlation can be used in the upscaling of measured relative permeability curves from laboratory scale to reservoir grid scale.

Inversion analysis of oil drainage test data
In this section an integrated approach to history matching of the oil drainage test data, which were carried out by unsteady-state method using glass bead models, is presented. The inversion algorithm includes a forward numerical modeling of fluid flow to obtain experimentally derived relative permeability curve. The history matching of the displacement process was optimized by using the Levenberg-Marquardt algorithm to minimize the error between the simulated and experimental production data (oil and water). In this algorithm, Corey-type power law is used to create relative permeability curves during the optimization procedures.
In the first step, a relative permeability curve is created by using Corey's form where the phase relative permeability is function of its saturation as described below: where k rw is the water relative permeability, k ro is the oil relative permeability, s wr is the residual water saturation, s or is the residual oil saturation, and a, b, c, and d are the parameters, which are adjusted to achieve an acceptable history match. An initial value is given to each of these controlling parameters (a, b, c, and d) to construct the first set of relative permeability curves. Notable to mention, the relative permeability end points of water and oil are determined by the experiments. Then, the simulation output (produced volumes of oil and water with time) is used to calculate the least square function, J, as follows: where Q obs i is the observed volume of produced fluid, Q sim i is the simulated volume of produced fluid, and N is the total number of recovery observations to be history matched.
In step two, the parameter ''a'' is modified by a magnitude of e (e = 0.001) to estimate the relative permeability curve which is then used to simulate two-phase fluid flow to calculate new fluid volumes (oil and water). These water and oil volumes are subsequently used in Eq. 10 to calculate a new least square function, J a , and then rJ a as follows: In step three, the sensitivity coefficient matrix, D, is calculated as: Now the first part of the left-hand side of the Levenberg-Marquardt algorithm is obtained and is known as Hessian matrix which is given by: In step four, similar procedures are repeated for other controlling parameters b, c, and d until D T D matrices are complete.
In step five, the full formulation of left-and right-hand sides of Levenberg-Marquardt algorithm is obtained. The algorithm is solved for the improvement term Dx k , and the equation used in this algorithm is as follows: where k is the stability factor, and Dx k is used to update the control parameters (a, b, c, and d) as: where a is the step size. Steps from one to five are repeated until Dx k be very small, and error between the observed and calculated volume of produced fluid is minimized.

Testing of inversion algorithm for two-phase flow
A 3D reservoir model with one fracture in the direction of fluid flow with dimensions of (10 m 9 10 m 9 10 m) is constructed to test the history matching optimization algorithm. A 3D finite element mesh is created, and the reservoir model is initialized by rock and fluid properties. The permeability of matrix and fracture is set to 0.1 md and 1D, respectively, and water and oil viscosity are set to 1 cp and 0.7 cp, respectively. A set of synthetic relative permeability curve is used as a target set. The created 3D model is saturated by water. Oil phase is injected from one end of the model at constant flow rate to displace the water phase, and the produced water volume (target volume) is calculated against time using the multiphase numerical model. The relative permeability curves for oil and water phases are changed, and then, new simulation outputs are obtained (simulated volume). The optimization history matching algorithm is used in this step to iterate the relative permeability until the error between the target and simulated volume of produced water is minimized. Figure 7 shows the three relative permeability curves include: target curve, initial curve, and optimized curve. The target and final relative permeability curves show a good match confirming that the optimization algorithm is capable of predicting relative permeability curves. Figure 8 shows a comparison between the optimized and target volume (synthetic) of water produced with an error less than 2 %. As can be seen in Fig. 8, the oil production rate is fluctuated and this can be explained by the following: At the beginning of the production process, the water produced from the fracture and then the matrix starts to feed the fracture with water. Therefore, during this process the water production rate will be fluctuated and it depends on how quick the matrix feed the fracture with fluid.

Numerical simulation of two-phase flow in glass bead fracture model
The first step is to numerically simulate Fahad (2013) glass bead experimental model and then predict oil drainage test data (produced volume of water). Next is to compare the simulated data with that obtained from the glass bead model in the laboratory. Finally, the inversion algorithm is used to match the recorded volume of produced water to predict relative permeability curves. 3D mesh for the glass bead experimental models with single and multiple fractures is presented in Fig. 9a-c. The 3D mesh is generated for the same physical dimensions as that of the laboratory model (20 cm 9 10 cm 9 0.2 cm). Matrix is represented by tetrahedral element, while the fractures by triangle elements. Matrix permeability is set to 3.4D, whereas fracture permeability is set to 1 9 10 4 D (from Fahad 2013). Oil viscosity is 1.4 cp and density is 0.75 gm/cm 3 , while water viscosity is 1.002 cp and density is 0.998 gm/cm 3 . Oil is injected at one end of the model at constant rate of 2 cm 3 /min, and the produced fluids are collected at the other end of the model (see Fig. 10). The simulation runs are performed without considering the effect of capillary pressure as no capillary pressure curve is available for the glass bead model.

Results
In Figs It is also observed from the displacement profile that the displacement process is delayed in the numerical simulation than that in laboratory experiment (see Figs. 14,15,16). This effect might be caused by ignoring the capillary forces in numerical modeling. Another possible reason is that the estimated permeability (by cubic law) for the fracture is lower than the actual permeability of the glass bead laboratory model.

Discussion
The relative permeability curves for oil and water obtained with minimum error (optimized) between the simulated and experimentally obtained produced water volume at 2 cm 3 /min constant injection rate are presented in Figs. 17 and 18, respectively. These curves are for single fracture at 0°and 45°and 2 intersected fractures at 45°orientation angle. As can be seen from Figs. 17 and 18, in general for single fracture oil relative permeability is lower than that of multiple fractures, while water relative permeability for single fracture is higher than that of multiple fractures. This can be explained by the fact that a large volume of water passes through the multiple fractures during the displacement process making the relative permeability to water for multiple fractures lower. From the results it can be concluded that the numerical methodology developed as part of this paper is capable of simulating the experimental observation made by Fahad (2013). Fahad (2013) conducted a large number of experiments in the glass bead model, and based on the results of these experiments, the authors developed a correlation model as follows: where k ro is the oil relative permeability, k ref is the reference oil relative permeability which is k ro of the homogenous glass bead pack (Fahad 2013), n f is the number of fractures in the system, h is the fractures orientation toward the flow direction, and w is the angle between the fracture with a reference fracture (one fracture that exists in the system is assumed to be a reference).
Results of this correlation are compared with the relative permeability obtained by simulation studies and presented in Fig. 19. It can be seen from Fig. 19 that the simulated oil relative permeability and that obtained based on the derived correlation are in a good agreement. From the comparison it is clear that the correlation has the power to Volume of Water Produced cm 3

Time (min)
Target (Synthetic) Optimized Fig. 8 Comparison between target (synthetic) and optimized volume of the produced fluid by using target and optimized sets of relative permeability (presented in Fig. 7) obtain the relative permeability of fracture systems, but this correlation is limited to the glass bead size.

Conclusion
In this paper, we have presented a methodology of estimating the relative permeability curve from drainage tests data using our developed multiphase fluid flow simulation model and an inversion technique. The multiphase fluid flow model is formulated based on the three governing forces, gravity viscous, and capillary forces. The flow process in porous fractured rock is simulated using finite element technique under poroelastic frame work. The model is validated against first, the analytical solution by Buckley-Leverett, simulating the oil drainage test in a glass bead model by Fahad (2013). From the validation, it was observed that the simulation results are in a good agreement with that at each stage of validation. The integration between the developed multiphase fluid flow simulation model and the inversion technique has been implemented to estimate the relative permeability curve from the laboratory test data that has been carried out on a different glass bed models. The results show that the developed integrated methodology has the ability to estimate the relative   Fig. 13 Oil-water drainage profile for 29 intersected fracture system at 45°in the direction of fluid flow at 0.4 PVI permeability curves for fractured systems. Then, these results gave us a concept of how this developed integrated methodology will be used in the future for upscaling relative permeability curves of fractured reservoirs from laboratory scale to reservoir scale.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix: Derivation of multiphase flow equations in a poroelastic framework
The mathematical equations of multiphase fluid flow through matrix and fracture system are derived and then reformulated using finite element technique. The combining of mass and momentum balance equations is used in the derivation process with appropriate boundary conditions and a set of assumptions. These assumptions include: (1) The system is only for water and oil fluid flow, and (2) the fluid is isothermal and fluid phases are in a thermodynamic equilibrium state.

Mass conversation equation
The conversation equations of masses of fluid and solid for a representative element volume (V) in a porous medium will be derived separately. The mass of solid component of the medium can be expressed as a volume fraction of the porous solid as the following equation shows, where / is the sum of fracture and matrix porosities and q s is the solid density The mass conversation of solid constituent can be written as; Equation (18) can be simplified if the continuum mechanism is assumed to be revealed; Mass conversation for oil and water phases in rock and fracture systems can be derived similarly as follows: The mass of solid component of the medium is expressed as: The mass conversation of fluid constituent can be written as; where w stands for oil and water phases, S w is the saturation of fluid phase (oil and water), q w is the fluid density, and q w is the rate of fluid exchange between matrix and fracture system. Equation (21) can be written separately for oil and water fluid phases separately by assuming that the continuum mechanism is prevailing as follows: where U w and U o are the intrinsic velocities for water and oil phases in matrix and fractured system, respectively. The assumption of quasi-steady-state flux between the fracture network and matrix is used; therefore, flux rate can be expressed as: where k 1 is referring to the permeability of the matrix and p 1w ; p 2w are the pressure in the matrix and fracture system, respectively. Darcy velocities for oil and water phases in matrix and fracture system are defined as: From Eq. (25) intrinsic phase velocities can be expressed as: Substituting Eq. (26) into Eq. (23) and then the equation reformulated as follows: Expanding the derivatives of Eqs. (19) and (27) Equations (29) and (30) can be reformulated and written as follows: By considering the total derivative as: Equations (31) and (32) can be written by considering the total derivative form as: From Eq. (34), the expression of change of porosity with time can be written as: Substituting Eq. (36) into Eq. (35), the following equation can be got: Or.
The relationships between the change of fluid and rock densities with fluid and rock bulk modulus are given as: where K w and K o are the bulk modulus of water and oil, respectively. K s and K ns are the bulk modulus of solid rock. Substituting Eq. (39) into Eq. (38) gets the following equation: In order to simplify Eq. (40) and get the final form of the two-phase fluid flow equation, the following assumptions are used as: • The total derivative is considered as: Eq. (C.24) with Darcy's law: • Solid velocity is small and can be neglected where w is standing for oil and water phase, k rw is phase relative permeability, q w is phase fluid density, g is the gravitational acceleration, and h is the height above reference level.
The water-phase flow governing equation will be as follows: The oil-phase fluid flow governing equation can be written as: where P wm and P wf are the water pressures inside matrix and fractures, respectively.

Momentum balance equation
The relationship between total applied stresses r ij and intergranular (effective) stresses r 0 ij is given by: where a is pore pressure ratio factor, and d ij is the Kronecker delta.The linear constitutive relationships of the system can be expressed as: where D ijkl is the elasticity matrix. The equilibrium equation of motion for a solid can be defined as: where F is the vector of tractions applied on the body. The strain-displacement relationship is defined as: Rearranging Eq. 42 for water phase and introduce strain, then same procedure for oil phase: Finite element discretization In this section the finite element technique is used to derive the integral formulation of the derived coupled fluid flow and rock deformation equations through fractured system. In addition, the finite element technique is used to discretize the problem domain into nodes and elements. In this paper, four-node tetrahedral elements were used to represent the rock matrix in a 3D space, while the discrete fractures are represented by triangle elements in a 2D space. The value of the material properties is assumed to be constant within the element and allowed to vary from one element to the next. For the three-dimensional four-node tetrahedral element, the shape functions have the following form: where n; g; f ð Þare the element local coordinates system. In order to transform the element geometry from the global system coordinates (x, y, and z) to the local coordinates n; g; f ð Þ, the following equations were used as: where J is the Jacobin matrix and can be expressed as: Shape functions are used to obtain the variation of unknown variables within the element. These variables are approximated by using the interpolation function in finite element space as follows: where N is the corresponding shape function, u, P w , and P o are the nodal unknown variables. B is the strain displacement matrix and can be defined as: Equilibrium equation can be written in its general form as follows: where of is the load vector applied on the boundary and V is the volume of the element. Equation (57) can be written as: Average pore pressure (P) can be defined as: Substituting Eq. (58) into Eq. (59); Equation (60) can be modified by inserting the capillary pressure as: Discretization form for water-phase flow equation through fractured porous media is given below as: The process is repeated for oil-phase flow in fractured porous media. All of the previous equations are used for fluid flow through discrete fractures, but in 2D space the final flow equation through fracture network and matrix is given as: where X f represents the fracture part of the domain as a 2D entity, and X m represents matrix domain and X is the entire domain.