Geometric model of the fracture as a manifold immersed in porous media

TL;DR

This paper introduces a geometric model for fluid flow in fractured porous media, using Riemannian manifolds to accurately represent complex fracture geometries and develop reduced models for efficient large-scale simulations.

Contribution

It develops a novel geometric framework modeling fractures as manifolds immersed in space, coupling it with porous media flow, and validates reduced models for complex fracture geometries.

Findings

The reduced model closely approximates the original geometric model.

The approach effectively handles complex, variable-thickness fractures.

The model facilitates large-scale reservoir simulations with complex fracture networks.

Abstract

In this work, we analyze the flow filtration process of slightly compressible fluids in porous media containing man made fractures with complex geometries. We model the coupled fracture-porous media system where the linear Darcy flow is considered in porous media and the nonlinear Forchheimer equation is used inside the fracture. We develop a model to examine the flow inside fractures with complex geometries and variable thickness, on a Riemannian manifold. The fracture is represented as the normal variation of a surface immersed in $\mathbb{R}^3$. Using operators of Laplace Beltrami type and geometric identities, we model an equation that describes the flow in the fracture. A reduced model is obtained as a low dimensional BVP. We then couple the model with the porous media. Theoretical and numerical analysis have been performed to compare the solutions between the original geometric model and the reduced model in reservoirs containing fractures with complex geometries. We prove that the two solutions are close, and therefore, the reduced model can be effectively used in large scale simulators for long and thin fractures with complicated geometry.

Tables

Table 1. Table 1: Coupled domain 1: comparison of the diffusive capacities obtained from the original model and the Reduced Model I.

	0.01			0.05			0.1
	$P{I}_{PD}$	$P{I}_{R_1}$	error	$P{I}_{PD}$	$P{I}_{R_1}$	error	$P{I}_{PD}$	$P{I}_{R_1}$	error
0	0.035027478	0.03524709	6.27E-03	0.05875387	0.057609773	1.95E-02	0.076646162	0.071405654	6.84E-02
0.001	0.034974637	0.03519495	6.30E-03	0.058705992	0.057564145	1.95E-02	0.076612623	0.071375039	6.84E-02
1	0.027792767	0.028039564	8.88E-03	0.041956574	0.041597979	8.55E-03	0.058528659	0.055027505	5.98E-02
10	0.025385268	0.025591167	8.11E-03	0.030338957	0.03057647	7.83E-03	0.037144807	0.036224497	2.48E-02
50	0.024668646	0.024847989	7.27E-03	0.026591578	0.027038853	1.68E-02	0.02924881	0.029500624	8.61E-03
100	0.024491447	0.024661886	6.96E-03	0.025662443	0.026149571	1.90E-02	0.027264275	0.027825031	2.06E-02

Table 2. Table 2: Coupled domain 2: comparison of the diffusive capacities obtained from the original model and the Reduced Model I.

	0.01			0.05			0.1
	$P{I}_{PD}$	$P{I}_{R_1}$	error	$P{I}_{PD}$	$P{I}_{R_1}$	error	$P{I}_{PD}$	$P{I}_{R_1}$	error
0	0.149934761	0.150403058	3.12E-03	0.184953216	0.182806385	1.16E-02	0.209349275	0.201457857	3.77E-02
0.001	0.149925428	0.150393766	3.12E-03	0.184940273	0.182794318	1.16E-02	0.209339858	0.201449376	3.77E-02
1	0.145535944	0.145975567	3.02E-03	0.176042698	0.174430344	9.16E-03	0.201650306	0.194507842	3.54E-02
10	0.140394828	0.140694173	2.13E-03	0.157643454	0.156802125	5.34E-03	0.177091286	0.172336013	2.69E-02
50	0.138128196	0.138326124	1.43E-03	0.147357331	0.146777702	3.93E-03	0.158752511	0.155725614	1.91E-02
100	0.137503318	0.137668886	1.20E-03	0.144376243	0.143869846	3.51E-03	0.153026912	0.150547714	1.62E-02

Equations

$$\frac{\mu }{k}\mathbf{v}+\beta |\mathbf{v}|\mathbf{v}=-\nabla p,$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{v})=0.$$

$${\rho }^{\prime }={\gamma }^{-1}\rho \left(\rho ={\rho }_0{\exp }^{{\gamma }^{-1}(p-p_0)}\right),$$

$$\frac{\partial p}{\partial t}=-\gamma \nabla \cdot \mathbf{v},$$

$$-\nabla p-\frac{\mu }{k}\mathbf{v}-\beta |\mathbf{v}|\mathbf{v}=0.$$

$$\mathbf{v}={\mathbf{v}}_{\beta }=-f_{\beta }\left(\Vert \nabla p\Vert \right)\nabla p,$$

$$f_{\beta }\left(\Vert \nabla p\Vert \right)=\frac{2}{\alpha +\sqrt{{\alpha }^2+4\beta \Vert \nabla p\Vert }}.$$

$$\displaystyle\Big{(}f_{\beta}(\|\eta_{1}\|)\eta_{1}-$$

$$\displaystyle\geq\frac{1}{2}f_{\beta}\big{(}\max(\|{\eta}_{1}\|,\|{\eta}_{2}\|)\big{)}\|{\eta}_{1}-{\eta}_{2}\|^{2}\,.$$

$$\frac{\partial p}{\partial t}=\gamma \nabla \cdot \left(f_{\beta }(\Vert \nabla p\Vert )\nabla p\right).$$

$$J_p(t)=\frac{{\int }_{{\Gamma }_w}\mathbf{v}\cdot \mathbf{n}ds}{{\bar{p}}_{\Omega }(t)-{\bar{p}}_w(t)},$$

$${\gamma }^{-1}\frac{\partial p}{\partial t}=\nabla \cdot \left(f_{\beta }(\Vert \nabla p\Vert )\nabla p\right),$$

$${\int }_{{\Gamma }_w}{f}_{\beta }\left(\Vert \nabla p\Vert \right)\frac{\partial p}{\partial \mathbf{n}}ds=-Q,$$

$$\displaystyle\frac{\partial p}{\partial\mathbf{n}}\Big{\rvert}_{\Gamma_{out}}=0\,,$$

$$p(x,t_0)=p_0(x).$$

$$-\nabla \cdot \left(f_{\beta }(\Vert \nabla W\Vert )\nabla W\right)$$

$$\displaystyle W\big{\rvert}_{\Gamma_{w}}=0\,,$$

$$\displaystyle\frac{\partial W}{\partial\mathbf{n}}\Big{\rvert}_{\Gamma_{out}}=0\,.$$

$${\int }_{{\Gamma }_w}{f}_{\beta }(\Vert \nabla W\Vert )\frac{\partial W}{\partial \mathbf{n}}ds=-Q.$$

$$p(x,t)=W(x)-\gamma At+K,$$

$$-\nabla \cdot k_p\nabla W_p$$

$$-\nabla \cdot f_{\beta }(\Vert \nabla W_f\Vert )\nabla W_f$$

$${\mathbf{v}}_p\cdot {\mathbf{n}}_p$$

$$W_p$$

$${\mathbf{v}}_p\cdot {\mathbf{n}}_p$$

$$W$$

$$G=\left[\begin{array}{c}{G}_{ij}\end{array}\right]$$

$$\Vert \mathbf{\Phi }{\Vert }_{{}_G}=$$

$${\nabla }_{{}_G}\varphi$$

$$\Vert {\nabla }_{{}_G}\varphi {\Vert }_{{}_G}^2$$

$${\nabla }_{{}_G}\cdot \mathbf{\Phi }=\sum _i\frac{1}{\sqrt{\left|G\right|}}{\partial }_i(\sqrt{\left|G\right|}{\varphi }^i).$$

$${\Delta }_{{}_G}\varphi ={\nabla }_{{}_G}\cdot {\nabla }_{{}_G}\varphi =\frac{1}{\sqrt{\left|G\right|}}\sum _{i,j}{\partial }_i\left(G^{ij}\sqrt{\left|G\right|}{\partial }_j\varphi \right).$$

$$f_{\beta }\left(\Vert {\nabla }_{{}_G}p{\Vert }_{{}_G}\right)=\frac{2}{\alpha +\sqrt{{\alpha }^2+4\beta \Vert {\nabla }_{{}_G}p{\Vert }_{{}_G}}},$$

$$\mathbf{v}=-f_{\beta }\left(\Vert {\nabla }_{{}_G}p{\Vert }_{{}_G}\right){\nabla }_{{}_G}p,$$

$$\alpha \mathbf{v}+\beta \Vert \mathbf{v}{\Vert }_{{}_G}\mathbf{v}=-{\nabla }_{{}_G}p.$$

$$\left(\beta \sqrt{\left({\nabla }_{{}_G}p\right)G\left({\nabla }_{{}_G}p\right)}{f}_{\beta }^2+\alpha f_{\beta }-1\right){\nabla }_{{}_G}p=0.$$

$$\left(\beta \Vert {\nabla }_{{}_G}p{\Vert }_{{}_G}{f}_{\beta }^2+\alpha f_{\beta }-1\right){\nabla }_{{}_G}p=0.$$

Full text

We develop a model to examine the flow inside fractures with complex geometries and variable thickness, on a Riemannian manifold. The fracture is represented as the normal variation of a surface immersed in $\mathbb{R}^{3}$. Using operators of Laplace Beltrami type and geometric identities, we model an equation that describes the flow in the fracture. A reduced model is obtained as a low dimensional BVP. We then couple the model with the porous media.

Theoretical and numerical analysis have been performed to compare the solutions between the original geometric model and the reduced model in reservoirs containing fractures with complex geometries. We prove that the two solutions are close, and therefore, the reduced model can be effectively used in large scale simulators for long and thin fractures with complicated geometry.

1 Introduction

Fractured reservoir modeling is a multi-scale and multi-physics complex problem, where analysis and simulation of this complex system require a deep understanding of the physical processes that describe the flow coupling at different scales. We consider flows in the porous media which include fractures of the same dimension as the reservoir itself. In practice fracture domain is very “long” compared to its width, and require a special care from analytical and numerical point of view. Moreover, the geometry of the fracture domain is a significant characteristic. Effects of complex geometric features of fractures such as variable thickness (which is very small) and the curvature of the fracture’s boundary are usually not taken into consideration during the simulations, which prevents successful reservoir modeling and leads to errors in forecasting the reservoir performance.

The traditional approach in the reservoir modeling is to simplify the problem by treating the fracture as a 1-D sink with pressure on the fracture equals to the value of the pressure on the well or with given flux, and the flow occurs in the porous media surrounding the fracture. Mathematically, this is formulated as a Dirichlet or Neumann boundary condition on the fracture for pressure function, which satisfies an elliptic or parabolic equation of second order. This assumption leads to a significant overestimation of the fracture capacity.

In our previous paper [19], we consider the nonlinear nature of the flow inside the fracture and couple it with the flow in the porous media. The main goal of this paper is to incorporate complex geometric features of the fractures into the modeling and explore the pressure distribution of the flow. We develop a framework based on methods of differential geometry to model the fractures with complicated geometries. In our approach, we formulate the fracture as a 3-D manifold immersed in porous media and we introduce a reduced dimensional nonlinear flow equation inside the fracture.

We consider slightly compressible fluid flows in the inhomogeneous reservoir fracture domain. Due to the heterogeneity of the fractured porous media, the velocity of the flow inside the fracture is much higher than inside porous media. It was observed that, due to high velocity inside fractures, the inertial effects become significant. Therefore the relation between the flow velocity and the pressure gradient inside the fractures deviates from Darcy’s law [8, 14, 9].

We consider Darcy’s law to model the flow in porous media and a generalized nonlinear Forchheimer equation to describe the flow inside the fracture. It is evident that the fluid mostly flows towards the fracture first and then transports to the well along the fracture. Consequently, the total production of the hydrocarbon depends on the capacity of the fracture to take-in fluid from the reservoir. This capacity depends on the geometry of the fracture and its conductivity, and it is characterized by the productivity index of the reservoir fracture system.

In past we proved that an integral functional called diffusive capacity, defined as the total flux on the well surface divided by the pressure drawdown (difference between average pressures of the reservoir domain and the well boundary) is mathematically equal to productivity index [7]. We use the diffusive capacity to characterize the reservoir performance.

We formulate a model to investigate the flow inside fractures with complex geometries. In particular, the fracture is represented as a perturbation across the fracture thickness in the normal direction to the barycentric surface immersed in $\mathbb{R}^{3}$, with its naturally induced Riemannian metric. Moreover, we formulate the flow equation inside the fracture using the first and second fundamental forms of the surface, geometric identities for Gaussian curvature and mean curvature, and corresponding operators of the Laplace-Beltrami type. On the boundary of the domain of the flow (union of porous media and fracture) we impose mixed boundary conditions. On the well we impose Dirichlet condition. This coupled fluid flow is impossible to solve numerically, since the thickness of the fracture is of the $10^{6}$ smaller than the length of the fracture. Therefore we introduce a reduced model for the flow inside the fracture domain. Moreover, we obtain further simplified models with further assumptions (that can be utilized depending on the size of the fracture thickness and other physical factors).

We theoretically and numerically investigate the difference between the solutions of the actual model and the reduced models. We confirm the successful implementation of the models, by proving that the solutions of the reduced models are close to the solutions of the actual model, for realistic values of fracture thickness.

Controlling the shape of the fractures in geological reservoirs is challenging. Therefore, our method can be applied mostly for simple fracture geometries. However, the geometric method and the analysis we introduce in this paper are valuable tools in modeling micro fluidic flows and blood flows in arteries and veins [6, 23]. Moreover, we believe that this methodology can be served as a foundation for reservoir engineers to model fractures in the future.

2 Formulation of the problem and preliminary results

In this section, we summarize important preliminary results on Darcy-Forchheimer equations from our paper [19].

2.1 Reservoir modeling

Mathematical framework of the reservoir modeling is based on Darcy-Forchheimer equation, the continuity equation and the state equation [7, 9, 18]. Among various methods that demonstrate non-Darcy case, non-linear Forchheimer equation is widely utilized [13, 12, 14, 21].

The velocity vector field $\mathbf{v}$ and the pressure $p$ in porous media are related by the Forchheimer equation given by

[TABLE]

where $k$ is the permeability, $\mu$ is the viscosity and $\rho$ is the density of the fluid. This describes the momentum conservation of the flow.

Remark 1**.**

In our intended application the parameters of the reservoir and the fracture are isotropic and space dependent. Namely, $k=k_{p}$ and $\beta=0$ in the porous block, and $k=k_{f}$ and $\beta\neq 0$ in the fracture.

The continuity equation takes the form

[TABLE]

For slightly compressible fluids, the state equation is given by

[TABLE]

where $\gamma$ is the compressibility constant of the fluid.

In natural reservoirs, the dissipation in the porous media is dominant [18]. We assume that the permeability coefficient $k$ is very small. Moreover, for many slightly compressible fluids $\gamma$ is of order $10^{8}$. [2, 9, 8]

With these constraints, the system can be rewritten as [19]

[TABLE]

Darcy-Forchheimer equation: The velocity vector field $\mathbf{v}$ can be uniquely represented as a function of the pressure gradient as follows.

[TABLE]

where $\alpha=\frac{\mu}{k}.$

Eq.(6) is referred as Darcy-Forchheimer equation. The velocity defined in Darcy Forchheimer equation (6) with $f_{\beta}$ defined by Eq.(7), solves the Forchheimer equation (5) [19].

The coefficient $\beta$ in the non linear term of the Eq.(6) does not depend on pressure. Changes of density of slightly compressible fluids have a minor impact on changes in the coefficient $\beta$ [9].

Remark 2**.**

As $\beta\to 0$ the Darcy-Forchheimer equation reduces to Darcy equation.

Function $f_{\beta}$ has important monotonic properties.

Lemma 2.1**.**

For $f_{\beta}(\|\eta\|)$ defined as above, the function $f_{\beta}(\|\eta\|)\|\eta\|$ is strictly monotonic on bounded sets. More precisely,

[TABLE]

The proof follows from Lemma 2.4 in Ref. [17] with $a=1/2$.

Lemma 2.2**.**

Let $f_{\beta}$ be defined by the formula (6). With the above assumptions, the pressure function $p$ satisfies the quasi linear parabolic equation

[TABLE]

Substituting Eq. (6) to the Continuity equation (4), Eq. (9) can be obtained. For more details see [20]

2.2 Diffusive Capacity and Pseudo-Steady State regime (PSS)

System (4) - (5) characterizes the fluid flow of well exploitation in a reservoir. Analogous to the reservoir engineering concept of productivity index (PI) that is used to measure the capacity and the performance of the well, we introduce the mathematical notion called the diffusive capacity.

Let $\Omega$ be a bounded reservoir domain bounded by the exterior no flux boundary $\Gamma_{out}$ and the well surface $\Gamma_{w}$. Let $\mathbf{n}$ be the outward unit normal on the piece wise smooth surface $\Gamma_{w}$.

Definition 2.3**.**

Diffusive Capacity.

Let the pressure $p$ and the velocity $\mathbf{v}$ form the solution for the system (4) - (5) in $\Omega$, with impermeable boundary condition \mathbf{v}\cdot\mathbf{n}\,\big{\rvert}_{\Gamma_{out}}=0. Then the diffusive capacity is defined by

[TABLE]

where, $\bar{p}_{\Omega}(t)-\bar{p}_{w}(t)$ is called the pressure draw down (PDD) on the well, $\bar{p}_{\Omega}(t)=\frac{1}{\left|\Omega\right|}\int_{\Omega}p\,d\Omega,\,\bar{p}_{w}(t)=\frac{1}{\left|\Gamma_{w}\right|}\int_{\Gamma_{w}}p\,ds,$ $\left|\Omega\right|$ is the volume of the reservoir and $\left|\Gamma_{w}\right|$is the area of the well.

It has been observed on field data that when the well production rate $Q(t)=Q=const$, the well productivity index stabilizes to a constant value over time [22].

Definition 2.4**.**

PSS regime.

Let the well production rate Q be time independent: $\int_{\Gamma_{w}}\mathbf{v}\cdot\mathbf{n}\,ds=Q.$ The flow regime is called a pseudo - steady state (PSS) regime, if the pressure draw down (PDD) = $\bar{p}_{\Omega}(t)-\bar{p}_{w}(t)$ is constant.

Corollary 1**.**

For a PSS regime the diffusive capacity/PI is time invariant.

2.2.1 PSS solution for the Initial Boundary Value Problem.

Let $\Omega$ be a bounded reservoir domain with impermeable exterior boundary and $p(x,t_{0})=p(x_{0})$ be the given initial pressure in the reservoir. Assume that the well is operating under time independent constant rate of production Q and the initial reservoir pressure is known.Then the IBVP that models the oil filtration process can be formulated as

[TABLE]

Since the boundary condition on $\Gamma_{w}$ is a single integral condition for the total flux, the IBVP (11)-(14) has infinitely many solutions. So, we constrain the solution to an auxiliary problem by assuming the solution to be constant in space on $\Gamma_{w}$. A more general situation has been considered in [5].

2.2.2 Steady state auxiliary BVP

Let $Q$ be the rate of production. Assume that the boundary $\partial\Omega$ is smooth. Let $W$ be the solution of the auxiliary steady state BVP:

[TABLE]

Through integration by parts we can obtain

[TABLE]

Proposition 1**.**

Let $W(x)$ be the solution of the auxiliary problem (15) - (17), then

[TABLE]

where $A=\frac{Q}{\left|\Omega\right|}$, solves the IBVP (11) - (14). For this solution we have the Pseudo Steady State and therefore the diffusive capacity is constant.

Proof.

Proposition 1 is followed by substituting (19) to Eq.(11) and verifying boundary and initial conditions (12) - (14). ∎

The constant $K$ in Eq. (19) is a measure of the initial oil reserve. With Definition (7) and monotonicity property in Lemma 2.1, the BVP (15)-(17) has a unique weak solution $W(x)$ belonging to $\mathbb{W}^{1,\frac{3}{2}}(\Omega)$ [4]. However, hereafter we only consider solutions $W(x)$ that satisfy $W\in C^{2}(\bar{\Omega})$ and $W,\nabla W\in C^{2}(\bar{\Omega})$ [17].

2.3 Fractured reservoir modeling

In this section, we introduce a fracture to the reservoir domain. We model the oil filtration process with a constant rate of production $Q$, for a fractured reservoir system. (For more details please see [19].) Consider a fractured reservoir domain where the exterior boundary $\Gamma_{out}$ is impermeable. Let $\Omega_{p}$ be the porous media domain, $\Omega_{f}$ be the fracture domain, $\Gamma_{f}$ be the boundary between the fracture and the porous media, and $\Gamma_{f_{out}}$ be the extreme of the fracture.

Let $\mathbf{v}_{p},W_{p},k_{p}$ be the velocity, the pressure and the permeability of the flow in porous media, respectively. Let $\mathbf{v}_{f},W_{f}$ be the velocity and the pressure of flow inside the fracture, respectively. Let $\mathbf{n}_{p}$, and $\mathbf{n}_{f}$ be the unit outward normal vectors to the porous medium and the fracture, respectively.

The auxiliary BVP (15) - (17) for the above mentioned reservoir-fracture system can be modeled as:

[TABLE]

Eqs. (23) and (24) assure the continuity of the solutions and the continuity of the fluxes across the interface $\Gamma_{f}$, respectively. Notice that we consider the linear Darcy law inside the reservoir while the non-linear Forchheimer equation is considered inside the fracture.

3 Geometric modeling of the fracture

Fractures in porous media have very complicated geometry. The domain of the fracture is very long compared to its thickness. The fracture thickness is changing along the length and the fracture boundary has curvature. As discussed in the introduction, neglecting those geometric features lead to over estimation of the fracture capacity. In this section, we employ methods in differential geometry to model the fractures with complex geometries. Then we obtain an equation for flow inside these fractures.

First, we introduce some definitions [10, 15, 16] that we will use in the modeling and obtain some preliminary results.

Definition 3.1**.**

(See [10]) An immersion is a differentiable function between differentiable manifolds whose derivative is everywhere injective. Explicitly, $f:M\to N$ is an immersion if $D_{p}f:T_{p}M\to T_{f(p)}N$ is an injective function at every point $p$ of $M$ (where $T_{p}X$ denotes the tangent space of a manifold $X$ at a point $p$ in $X$).

Let $\mathbf{R}:D\subset\mathbb{R}^{3}\to\mathbb{R}^{3}$; $\mathbf{R}(u,v,w)=\left\langle X(u,v,w),Y(u,v,w),Z(u,v,w)\right\rangle$ represents an immersion of a three dimensional object $M$ in $\mathbb{R}^{3}$; $M=R(D)$ where $D$ is an open simply connected domain.

Definition 3.2**.**

The induced metric associated to $\mathbf{R}(u,v,w)$ is defined as

[TABLE]

where $\langle\cdot,\cdot\rangle$ represents the inner product in $\mathbb{R}^{3}$.

We will denote by $(M,G)$, the induced Riemannian manifold with metric $G$ on $M$. For the next definitions, please see [15].

Definition 3.3**.**

The norm of a vector $\mathbf{\bm{\Phi}}=\langle\phi^{1},\phi^{2},\phi^{3}\rangle$ on the manifold is defined by

[TABLE]

Definition 3.4**.**

Let $G^{-1}=G^{ij}$ be the inverse of G. The gradient of a differentiable function $\phi$ on the manifold is defined by

[TABLE]

Note:

In local coordinates, the $j^{th}$ component of the vector field $\nabla_{{}_{G}\,}\phi=\left(\nabla_{{}_{G}\,}\phi\right)^{j}=\sum_{i}G^{ij}{\partial_{i}}\phi\\ =\left(G^{-1}\nabla\phi\right)^{j}$ with corresponding basis $\left\{\partial_{j}\right\}_{j=\{1,2,3\}}=\left\{\frac{\partial}{\partial x_{j}}\right\}_{j=\{1,2,3\}}$

Using Eqs. (27) and  (28), we obtain the norm of the gradient on the manifold

[TABLE]

Definition 3.5**.**

Let $|G|$ be the determinant of G. The divergence of a vector field $\bm{\Phi}$ on the manifold is defined by

[TABLE]

Using Eqs. (28) and (30), the Laplace-Beltrami operator applied to a scalar $\phi$ is given by

[TABLE]

Next, we obtain the Darcy-Forchheimer equation on the manifold $M$.

Lemma 3.6**.**

Let

[TABLE]

where $p$ is the pressure. Then the velocity defined by the Darcy Forchheimer equation on the manifold M

[TABLE]

solves the Forchheimer equation on M ,

[TABLE]

Proof.

Using Eq. $\eqref{DF_M}$, Eq. $\eqref{F_M}$ can be rewritten as

[TABLE]

Using Eq. (29), we have $\sqrt{\left(\nabla_{{}_{G}\,}p\right)G\left(\nabla_{{}_{G}\,}p\right)}=\|\nabla_{{}_{G}\,}p\|_{{}_{G}\,}$ and therefore

[TABLE]

We need to show that (32) solves the above equation. For $\nabla_{{}_{G}\,}p=\bm{0}$ it is true for any $f_{\beta}$. For $\|\nabla_{{}_{G}\,}p\|_{{}_{G}\,}\neq 0$, then $f_{\beta}$ should satisfy

[TABLE]

The positive root is given by

[TABLE]

or equivalently

[TABLE]

which is obtained multiplying both denominator and numerator of(38) by

[TABLE]

∎

3.1 Geometric model of the fracture as a 3-D manifold

In this section, we model the fracture as a 3-D manifold immersed in porous media. In our approach, we formulate the fracture as a parametrized Riemannian hypersurface. (Please see [20] for more details.) We describe the fracture as the normal variation of a surface $\mathbf{r}(u,v)$ immersed in $\mathbb{R}^{3}$ given by,

[TABLE]

where

[TABLE]

$\mathbf{n}$ is the outward unit normal vector to the surface such that

[TABLE]

and $\lambda\in[-h(u,v),h(u,v)]$. Here $2h(u,v)$ represents the thickness of the fracture.

Let $g=g_{ij}$ be the first fundamental form of $\mathbf{r}(u,v)$ given by

[TABLE]

and second fundamental form of $\mathbf{r}(u,v)$ is given by

[TABLE]

Since

[TABLE]

for ${i,j}=\{u,v\}$, the coefficients of the first fundamental form can be rewritten as

[TABLE]

[TABLE]

Let $|g|$ be the determinant of $g$, $K$ be the Gaussian curvature of $\mathbf{r}(u,v)$, and $H$ be the mean curvature of $\mathbf{r}(u,v)$ given by [10]

[TABLE]

Let $|G|$ be the determinant and $G^{-1}=G^{ij}$ be the inverse matrix of $G$, respectively. We have

[TABLE]

With this metric $G$ and the coefficients defined above, next, we obtain the equation for pressure of the flow inside the fracture domain.

3.2 Flow equation inside the fracture

Now we model the pressure distribution of a nonlinear flow inside the fracture domain $\Omega_{f},$ which is defined as a general manifold described above. The boundary of the domain of the flow is split as $\partial{\Omega_{f}}=\Gamma_{w}\cup\Gamma_{f}^{\pm}\cup\Gamma_{out}^{\pm}\cup\Gamma_{f_{out}}.$

[TABLE]

Schematically the domain of the flow with its boundaries is presented in the Figure 2.

Using Eqs. (21), (29) and (31) we have the pressure of the flow inside $\Omega_{f}$ given by the equation

[TABLE]

or equivalently

[TABLE]

where

[TABLE]

$i,j=\{1,2,3\}\,\,\text{and}\,\,\partial_{1}=\partial_{u},\partial_{2}=\partial_{v}\,,\,\partial_{3}=\partial_{\lambda}.$ Let $\mathbf{n}_{f}^{\pm}$ be the unit outward normal vector to the top/bottom boundary $\Gamma_{f}^{\pm}$ and $\mathbf{n}$ be the unit outward normal vector to the boundaries $\Gamma_{f_{out}}$ and $\Gamma_{out}^{\pm}$. Following mixed boundary conditions are imposed on the boundaries: Dirichlet condition on the well $\Gamma_{w}$, flux conditions on the top and the bottom boundaries $\Gamma_{f}^{\pm}$, no flux condition on outer boundaries $\Gamma_{f_{out}}$ and $\Gamma_{out}^{\pm}.$ Moreover, for simplicity we assume no flow in $v$ direction. Namely,

[TABLE]

3.3 Reduced model of the flow inside the fracture

In this section we simplify the original flow equation inside the fracture (given by Eq. (47)) and obtain a reduced dimensional model for the flow in the fracture.

Using the fact that $G^{23}=G^{32}=0,\,G^{33}=1$ and integrating Eq. (47) over the thickness of the fracture, we get

[TABLE]

where $i,j=\{1,2\}$.

Proposition 2**.**

Assume that the boundary conditions (48)-(50) hold. Then Eq. (51) can be reduced to

[TABLE]

where $i,j=\{1,2\}$,

[TABLE]

Proof.

Using the Leibniz rule, the first integral of the right hand side of Eq. (51) can be rewritten as

[TABLE]

where $h_{1}=h_{u}$ and $h_{2}=h_{v}$, and the second integral can be rewritten as

[TABLE]

Due to Eqs.  (28), (33) and (42)

[TABLE]

[TABLE]

Therefore, on $\Gamma^{+}_{f}$ we have,

[TABLE]

Similarly we can show that on $\Gamma^{-}_{f}$,

[TABLE]

Hence,

[TABLE]

Therefore, using Eqs.(59), (53) and (54), Eq. (51) can be rewritten as,

[TABLE]

∎

The thickness of the fracture is several orders smaller compared to the length of the fracture. Therefore it is reasonable to assume that the flow inside the fracture in the normal direction to the barycentric surface is negligible. With that assumption, we obtain the reduced model for the flow pressure in the fracture.

Proposition 3**.**

Reduced Model I. Let all the conditions of proposition 2 are satisfied. Assume that the gradient of $W$ is independent of $\lambda.$ Then the equation for pressure of the flow inside the fracture can be given by,

[TABLE]

where

[TABLE]

$i,j=\{1,2\}$**

We use numerical integration (3-point Gauss Quadrature Rule) to evaluate integrals $L^{ij}$ for $i,j\in\{1,2\}$ and $A$.

3.4 Reduced model of the fracture: when the metric $G$ does not depend on $\lambda$

For fractures with small thicknesses it is reasonable to assume that the solution of the flow equation in the direction of the thickness does not change, namely, we assume that the solution $W$ does not depend on the parameter $\lambda$ and consequently the first fundamental form of $\mathbf{R}(u,v,\lambda)$, $G$, does not depend on $\lambda$ as well. Therefore, next we further simplify the Reduced Model I.

Proposition 4**.**

Reduced Model II. Let all the conditions of theorem 3 are satisfied. In addition assume that the first fundamental form of $\mathbf{R}(u,v,\lambda)$ depends on the function $\mathbf{r}(u,v)$ only (does not depend on $\lambda$):

[TABLE]

Then the pressure inside the fracture is subjected to the equation:

[TABLE]

$i,j=\{1,2\}\,,$* where,*

[TABLE]

Proof follows from theorem 3.

Remark 3**.**

The reason why we consider two reduced models is the following. Reduced Model I is more comprehensive. It originates from the actual model under the assumption that the gradient of the pressure function $W$ does not depend on the parameter $\lambda$ which physically means that the velocity inside the fracture in the orthogonal direction to the barycentric surface is negligible. In this case, the coefficients of the equation $L^{ij}$ implicitly depend on $\lambda$ through integration.

Reduced Model II is a simplified version of the Reduced Model I, in which coefficients of the equation do not depend on $\lambda$. Consequently, the corresponding solution does not depend on $\lambda$. It is clear that, as the thickness of the fracture becomes big enough, the solution obtained from the Reduced Model II will significantly deviate from the actual model. In section 5, we investigate this numerically in detail.

4 Estimates for the difference between the solutions of the original model and the reduced models

In this section, we provide estimates for the difference between the solution of the original model and the solution of the Reduced Model II, when the fracture domain $\Omega_{f}$ is considered to be a foliation of a cylindrical surface in $\mathbb{R}^{3}$, namely

[TABLE]

We first investigate the difference between the solutions only inside the fracture with given fluxes, Theorem 4.5, and then we investigate the difference between the solutions in the coupled domain, Theorem 4.6.

Proposition 5**.**

Let all the conditions of proposition 4. Then the equation for pressure inside the fracture as a foliation of a cylindrical surface in $\mathbb{R}^{3}$is given by

[TABLE]

where

[TABLE]

Proof.

For a fracture with $z(u,v)=f(u)$ and thickness $2h(u)$, the barycentric surface of the fracture is given by $\mathbf{r}(u,v)=\left\langle u,v,f(u)\right\rangle$ and the inverse of the induced metric associated to $\mathbf{R}(u,v,\lambda)$ is given by

[TABLE]

Then the result follows from proposition 4. ∎

The analysis will be based on the following results.

Theorem 4.1**.**

*(See Ref.[11])

Let $M$ be a compact oriented Riemannian manifold of dimension $n$ with boundary $\partial M$. Then for all the vector fields $X$ and smooth functions $f$, the integration by parts formula on the manifold is given by*

[TABLE]

Definition 4.2**.**

On the manifold $M$, $L^{p}$ norm is defined as the following .

[TABLE]

Lemma 4.3**.**

For $f_{\beta}(\|\nabla_{g}W\|_{g})$ defined by Eq. (7), $1\leq q<2$,

[TABLE]

The proof of the above lemma can be obtained using the same arguments as for the case of $\mathbb{R}^{n}$ in Lemma III.11 in Ref [3] with $a=1/2$.

Lemma 4.4**.**

There exists a constant $C$ depending on $\Omega$, $Q$, $q^{+}(u)$ and $q^{-}(u)$ such that the corresponding basic profiles $W$ and $\bar{W}$ satisfy

[TABLE]

The proof of the above lemma can be obtained using the same arguments as for the case of $\mathbb{R}^{n}$ in Theorem V.4 in Ref [3] with $a=1/2$.

4.1 Estimates for the difference between the solutions inside the fracture $\Omega_{f}$.

We investigate the difference between the solutions of the actual model (i) and the reduced model (ii) on the 2-D manifold as defined below. Let $W$ be the solution of the actual model and $\bar{W}$ be the solution of the reduced model. Here $\nabla_{g}=<\frac{1}{1+f_{u}^{2}}\partial_{u},\partial_{\lambda}>$ and $\nabla_{g_{u}}=\frac{1}{1+f_{u}^{2}}\partial_{u}.$

1. (i)

Actual model: The flow equation is given by

[TABLE]

with the boundary conditions

[TABLE] 2. (ii)

Since the solution of the reduced problem is $\lambda$ independent, for comparison of the actual problem and the reduced one, we state the 1-D reduced problem as a 2-D one in the same domain $\Omega_{f}$.

Reduced model: The flow equation is given by

[TABLE]

with boundary conditions

[TABLE]

Theorem 4.5**.**

Let $W$ and $\bar{W}$ be the solutions of B.V.P.s (i) and (ii) respectively. Then,

[TABLE]

for some constant $C$.

Proof.

Subtracting Eq. (72) from Eq. (68), multiplying by $z=z(u,\lambda)=W-\bar{W}$ and integrating over the volume of the fracture we obtain

[TABLE]

Using Green’s formula on Riemannian manifolds (See Ref [11]) and due to the boundary conditions (70), (71), (74) and (75), the left hand side of the above equation can be rewritten as

[TABLE]

In the previous equations, as well as in the rest of the paragraph, $d\mu$ is the Riemannian volume element of the manifold $M$ and $d\tilde{\mu}$ is the area element of the boundary of the manifold $M$. Eq. (77) is equivalent to

[TABLE]

where,

[TABLE]

Consider $I_{1}$. By lemma 4.4 and 4.3 with $q=3/2$, there exist positive constant $C_{0}$ such that,

[TABLE]

Now consider $I_{3}$. Due to the boundary conditions (69) and since $\bar{W}$ is $\lambda$ independent, we have

[TABLE]

It then follows that

[TABLE]

Using Hölder and Cauchy inequalities we obtain

[TABLE]

Using Poincaré inequality, we have

[TABLE]

and

[TABLE]

Therefore,

[TABLE]

Combining Eqs. (82), (84) and (88), choosing $\varepsilon=\frac{C_{0}}{4}\left(\frac{3}{2}\right)^{\frac{4}{3}}$ and setting $C_{1}=\frac{1}{C_{0}}\left(\frac{2}{3}\right)^{\frac{4}{3}}$ yields

[TABLE]

Therefore, we have

[TABLE]

with $C=2C_{1}/C_{0}$. ∎

Remark 4**.**

From the theorem above, it follows that for a given fracture with thickness $h$, the difference between the solutions of the two problems can be controlled by the boundary data.

However, it should be noted that in the reservoir-fracture system as $h$ goes to zero, the fracture vanishes, and the oil flows mostly towards the well. Then as $h$ becomes smaller, $q^{+}$ and $q^{-}$ gets smaller as well, and therefore the individual velocities remain bounded.

Next, we investigate the coupled fractured-porous media domain. We show a much stronger result, under the condition that as $h\rightarrow 0$ the fluxes on the fracture boundary vanish with the same speed.

4.2 Estimates for the difference between the solutions in coupled fracture-porous media domain with linear isotropic flows.

Next, we provide estimates for the difference between the solutions in coupled domain. We consider half of a symmetric idealized fracture-reservoir domain depicted below and consider the flow to be linear isotropic inside the fracture. Let $\Omega_{p}$ and $\Omega_{f}$ be the porous media region and the fracture, respectively, with $\Omega=\Omega_{p}\cup\Omega_{f}$. Let $\Gamma_{f}$ be the top boundary of the fracture described by $\mathbf{R}(u,v,h)$, $\Gamma_{w}$ be the well boundary, $\Gamma_{out}$ be the outer boundary of $\Omega$, $\Gamma_{f_{out}}$ be the right extremum of the fracture. Let $\mathbf{n}$ and $\mathbf{n}_{f}$ be the outward unit normal on $\Gamma_{out}$ and $\Gamma_{f}$ respectively. We build the domain $\Omega_{f}$ such that $\mathbf{R}(u,v,\lambda)=\mathbf{r}(u,v)+\lambda\,\mathbf{n}_{f}$, with $0\leq\lambda\leq h$. Namely, we get the profile $\mathbf{r}(u,v)$ imposing $\mathbf{r}(u,v)+h\,\mathbf{n}_{f}=\mathbf{R}(u,v,h)$, with $\mathbf{R}(u,v,h)$ given.

Let $k_{p}$ and $k_{f}$ be the permeability of the porous media and the fracture respectively. Let $W_{i}$ be the flow pressure in the original problem and $\bar{W}_{i}$ be the flow pressure in the reduced problem with $i\in\{p,f\}$. Here $p$ denotes the porous media and $f$ denotes the fracture. Let $\bar{q}(u)$ be the flux coming into the fracture from the reservoir.

We investigate the difference between the solutions of the two problems defined below with $\nabla_{g}=<\frac{1}{1+f_{u}^{2}}\partial_{u},\partial_{\lambda}>$ and $\nabla_{g_{u}}=\frac{1}{1+f_{u}^{2}}\partial_{u}.$

1. (I)

The flow equations for the original problem are given by

[TABLE]

with boundary conditions

[TABLE] 2. (II)

The flow equations for the Reduced Model is given by

[TABLE]

with boundary conditions

[TABLE]

For simplicity, we drop the subscripts in $W_{p}$, $W_{f}$, $\bar{W}_{p}$ and $\bar{W}_{f}$ and use notations $W$ and $\bar{W}$ for the solutions in the original problem (I) and reduced problem (II) respectively. Each of the solution corresponds to their domain of integration.

We assume that, for a manifold with appropriate conditions on Riemannian metric, the following conjecture to be true.

Conjecture 1**.**

For $\bar{q}(u)$ given by Eq. (102), $\left|\dfrac{\bar{q}(u)}{h}\right|\leq C_{0}$ for some constant $C_{0}$ and for all $(u,v,\lambda)\in\Gamma_{f}$.

Henceforth $d\mu$ and $d\Omega$ are the elementary volumes of the manifold and the Euclidean domain respectively. $d\tilde{\mu}$ and $dS$ are the elementary areas of the manifold and the Euclidean domain respectively. By construction, the fracture domain $\Omega_{f}$ has the following property: the element of volume and the element of area of the fracture in Euclidean domain and on manifold are exactly the same, namely, $d\Omega=d\mu$ and $dS=d\tilde{\mu}$.

Theorem 4.6**.**

Let $W$ and $\bar{W}$ be the solutions of B.V.P.s (I) and (I) respectively. Then, under the assumption in conjecture (1) the following estimate holds for some constant $C$ that depends on $C_{0}$.

[TABLE]

Proof.

Let $z=z(u,v,\lambda)=W-\bar{W}$. Subtracting Eq. (98) and  (99) from Eq. (91) and  (92), respectively, we obtain

[TABLE]

and

[TABLE]

Then, after adding Eqs. (106) and (107), multiplying by $z$ and integrating over the volume of the domain we obtain,

[TABLE]

Using Green’s formula on Riemannian manifolds (See Ref [11]) and boundary conditions (95), (96), (97), (101), (103) and (104), the first and the second integrals in the left hand side of the above equation can be rewritten as

[TABLE]

and

[TABLE]

respectively.

Combining Eqs. (109), (110), boundary condition (94) and since \Big{<}f_{\beta}(\|\nabla_{g}\bar{W}\|_{g})\nabla_{g}\bar{W},\mathbf{n}_{f}\Big{>}_{g}=0, Eq. (108) can be rewritten as

[TABLE]

Due to boundary condition (102) and since $dS=d\tilde{\mu}$, we have

[TABLE]

Denote

[TABLE]

Then,

[TABLE]

Consider $I_{2}$. By Lemma 4.4 and 4.3 with $q=\frac{3}{2}$, there exist positive constant $C_{1}$ such that,

[TABLE]

Also, from conjecture (1), we have

[TABLE]

Using Hölder and Cauchy inequalities we obtain

[TABLE]

where $\left|\Omega_{f}\right|$ is the volume of the fracture. Then, using Poincaré inequality, we have

[TABLE]

Therefore,

[TABLE]

Combining Eqs. (114), (115), (119) we obtain

[TABLE]

Choose $\epsilon=\frac{1}{h^{2}}\frac{C_{1}}{C_{0}}\left(\frac{3}{2}\right)^{\frac{4}{3}}$. Then

[TABLE]

with $\tilde{C}=\dfrac{C_{0}^{2}}{4C_{1}}\left(\dfrac{2}{3}\right)^{\frac{4}{3}}$. Therefore,

[TABLE]

where $\left|\Gamma_{f}\right|$ is the surface area of the boundary $\Gamma_{f}$ and $C=\dfrac{\tilde{C}}{k_{p}}\left|\Gamma_{f}\right|^{\frac{2}{3}}.$ ∎

Remark 5**.**

From the theorem above, it can be observed that the estimate for the difference between the solutions of the original problem and the reduced problem goes to zero as $h$ going to zero, representing a much stronger estimate.

5 Numerical analysis and simulations

In this section, we present numerical results to demonstrate the validity of our approach. Namely, we show that the solutions of the Reduced Model I (3) and the Reduced Model II (4) are close to the solution of the original problem (46) for different fracture-reservoir geometries.

All the simulations have been performed using COMSOL Multiphysics software [1]. The grid size has been refined until changes in the pressure distribution between the previous and the next steps are negligible. The length and the thickness of the fracture is selected in relative units, and are dimensionless. Hydrodynamic parameters such as permeability and Forchheimer coefficient are numerically chosen without bonding to actual data of the porous media properties. In all simulations, the following parameters have been fixed: permeability in the porous media $k_{p}=0.01$, permeability inside the fracture $k_{f}=1$ and production rate $Q=1$.

5.1 Pressure distribution of the flow inside the fracture

In this section, we compare the pressure distributions of the flow obtained from the original model, the Reduced Model I and the Reduced Model II, inside the domain of fracture only, for different fracture geometries. For both Examples 5.1 and 5.2 we perform the following numerical simulations: First, we obtain the solution of the original model (46) inside the fracture, imposing zero Dirichlet boundary condition on the well, given flux boundary conditions on the top and bottom of the fracture, and zero Neumann boundary condition on the right end (namely, system (48)-(50) with $q^{\pm}=10$). Then, we solve the Reduced model I (61) and Reduced Model II (63) on the barycentric line of the fracture cross section, with zero Dirichlet boundary condition on the well and zero Neumann boundary condition on the right end.

Example 5.1**.**

We consider the fracture geometry with barycentric surface given by

[TABLE]

and constant thickness $2h(u)$. Since the solution of the problem does not depend on $v$, we solve our equations only on the cross section given in Figure 5.

In Figure 6, we compare the pressure distributions obtained from the original model (on the barycentric line of the fracture cross section) and the Reduced Models for $\beta=0.1$, $2h(u)=0.025$ (left) and $2h(u)=1$ (right). It is evident that the solutions are almost identical to each other when the thickness of the fracture is relatively small. When the thickness is large, it can be observed that the solution obtained from the Reduced Model I is very close to the solution of the original model, but the solution of the Reduced Model II deviates from the original one. This can be explained by the presence versus the absence of the parameter $\lambda$ in the Reduced Model I and II, respectively, as explained in Remark 3.

Example 5.2**.**

Next, we consider the fracture geometry with barycentric surface given by

[TABLE]

and variable thickness $2h(u)=0.2(2+0.5\sin(7u))$. Again, since the solution of the problem does not depend on $v$, we solve our equations only on the cross section given in Figure 7.

In this example, the thickness of the fracture changes depending on the parameter $u$. In Figure 8, we compare the pressure distributions obtained from the original model (on the barycentric line of the fracture cross section), the Reduced Model I and the Reduced Model II, for $\beta=0$ (Darcy) and $\beta=0.1$ (Forchheimer). The solutions of the Reduced Model I remain close to the solutions of the original model. However, the solutions of the Reduced Model II deviate from the solutions of the original model, especially in the Forchheimer case.

5.2 Diffusive capacity in the coupled fractured porous media domain

In this section, we calculate the diffusive capacities in the coupled fracture reservoir domain, using the original model (47) and the Reduced Model I (61).

The productivity index, which characterizes the well capacity to take-in hydrocarbons from reservoir, depends on the geometry of the fracture and its conductivity, and it is evaluated using the diffusive capacity of the well-reservoir-fracture system. Denote the diffusive capacities of the original model and the Reduced Model I by $PI_{PD}$ and $PI_{R_{1}}$, respectively. We compare the diffusive capacities obtained from the two models, as the amplitude of the thickness $H$ and Forchheimer coefficient $\beta$ changes, for different geometries of the fracture in the reservoir. Also, we evaluate the relative error as $\mbox{error}=\left|\frac{PI_{PD}-PI_{R_{1}}}{PI_{PD}}\right|.$ The goal here is to show that, in the fully coupled domain, the diffusive capacity calculated using the Reduced Model I is very close to the diffusive capacity calculated using the original model.

Example 5.3**.**

In this example, we numerically investigate the pressure distribution and the diffusive capacity of an infinite long reservoir, whose cross section is the rectangle $[-10,20]\times[-10,10]$. The well is modeled as an infinite long cylindrical surface with square cross section of side length $0.5$, centered on the y-axis and rotated by 60 degrees around it. The fracture barycentric surface is given by

[TABLE]

and variable thickness $2h(u)=H(2+0.5\sin(7u))$. The vector $\langle x_{c},y_{c},z_{c}\rangle=\langle\frac{1}{8},0,\frac{\sqrt{3}}{8}\rangle$ has been chosen so that the fracture starts from the center of the top-right face of the well. Since the solution of the problem does not depend on $y$, we solve our equations only on the cross section of the domain given in Figures 9 (left and right).

First, we couple the original flow equation inside the fracture with the flow in the porous media, by imposing the continuity of the solutions and the continuity of the fluxes across the fracture boundaries. Zero Dirichlet boundary conditions are imposed on the well. Zero flux boundary conditions are imposed on all the outer boundaries of the reservoir and on the right end of the fracture. Then, we solve the coupled system using the flow equations in the Reduced Model I with zero Dirichlet boundary condition on the well and zero Neumann boundary condition on the right end of the fracture.

Figure 9, presents the pressure distributions in the coupled domain obtained using the original model (left) and the Reduced Model I (right), for $H=0.1$ and $\beta=0.001$. The colors indicate that the fluid first converges towards the fracture and then flows towards the well. In Figure 10, we compare the pressure distributions obtained using the original model (on the barycentric line of the fracture cross section) and the Reduced Model I, for $\beta=0$ (Darcy) and $\beta=100$ (Forchheimer). It is evident that the solutions of the two models are very close to each other.

Table 1 presents the diffusive capacities in the coupled domain, obtained from the original model and the Reduced Model I, as $H$ and $\beta$ change. Moreover, the relative error, $\left|\frac{PI_{PD}-PI_{R_{1}}}{PI_{PD}}\right|$, has been reported.

Example 5.4**.**

Next, we evaluate the diffusive capacity in an infinite long reservoir whose cross section is the rectangle $[-15,15]\times[-10,10]$. The well is modeled as an infinite long cylindrical surface centered on the $y$-axis, with rectangular cross section of height $0.5$ and width $20$. The reservoir contains three fractures whose geometries are identical, and are described by the barycentric surface given by

[TABLE]

and variable thickness $2h(u)=H(2+0.5\sin(7u))$. In here

[TABLE]

and $A_{i}$ are rotation matrices given by

[TABLE]

Again, since the solution of the problem does not depend on $y$, we solve our equations only on the cross section of the domain given in Figures 11 (left and right).

First, we couple the original flow equation inside the fracture with the flow in the porous media, by imposing the continuity of the solutions and the continuity of the fluxes across the fracture boundaries. Zero Dirichlet boundary conditions are imposed on the well. Zero flux boundary conditions are imposed on all the outer boundaries of the reservoir and on the ends of the fractures. Then, we solve the coupled system using the flow equations in the Reduced Model I, with zero Dirichlet boundary condition on the well and zero Neumann boundary condition on the ends of the fractures.

Figure 11 presents the pressure distributions in the coupled fracture porous media domain, obtained from original model (left) and the Reduced Mode I (right), respectively, for $H=0.1$ and $\beta=1$. The colors indicate that the fluid first converges towards the fracture and then flows towards the well. Table 2 presents the diffusive capacities in the coupled domain, obtained using the original model and the Reduced Model I, as $H$ and $\beta$ change. Moreover, the relative error has been reported.

Both Examples 5.3 and 5.4 show similar observations. Tables 1 and 2 confirm that the diffusive capacities obtained from both models increase with increasing fracture thickness, for all $\beta$. At the same time, their values decrease, when the Forchheimer coefficient $\beta$ increases. This numerical results have a clear physical interpretation. Moreover, we can see the error is small for all values of $H$ and $\beta$. The error increases with $H$, while it has a slight dependence on $\beta$. Obtained results show that, as $H\approx 10^{-2}$, the errors are very small, and therefore the Reduced Model I can be effectively used in large scale simulators for long and thin fractures with complicated geometries.

6 Conclusion

In this paper, we investigated the flow filtration process of slightly compressible fluids in porous media containing fractures with complex geometries. We modeled the coupled fractured porous media system where the linear Darcy flow is considered in porous media and the nonlinear Forchheimer equation is used inside the fracture. Using methods in differential geometry we formulated the fracture as a manifold immersed in Porous media. The equation for pressure of the flow inside the fracture was modeled and then a reduced model, where the fracture is presented as a boundary inside porous media, was obtained. Theoretical and numerical results were obtained to prove the closeness of the solutions of the actual model and the solutions of the reduced model, both inside the fracture and in the coupled domain.

The main conclusion can be formulated as: for actual field data, where the thickness of the fracture is small compared to the length of the fracture, the reduced model can be effectively used in large scale simulators for long and thin fractures with complicated geometry.

Controlling the shape of the fractures in geological reservoirs is a challenging problem. Therefore our method can be applied mostly for simple fracture geometries. However, the geometric method and the analysis we introduced in this paper are valuable tools in modeling micro fluidic flows and blood flows in arteries and veins. Moreover, we believe that this methodology can be served as a foundation for reservoir engineers to model fractures in the future.

Acknowledgement

This work was supported by the National Science Foundation grant NSF-DMS 1412796.