Dimensional reduction of a fractured medium for a polymer EOR model

Dimensional reduction strategy is an effective approach to derive reliable conceptual models to describe flow in fractured porous media. The fracture aperture is several orders of magnitude smaller than the characteristic size (e.g., the length of the fracture) of the physical problem. We identify the aperture to length ratio as the small parameter 𝜖 with the fracture permeability scaled as an exponent of 𝜖. We consider a non-Newtonian fluid described by the Carreau model type where the viscosity is dependent on the fluid velocity. Using formal asymptotic approach, we derive a catalogue of reduced models at the vanishing limit of 𝜖. Our derivation provides new models in a hybrid-dimensional setting as well as models which exhibit two-scale behaviour. Several numerical examples confirm the theoretical derivations of the upscaled models. Moreover, we have also studied the sensitivity of the upscaled models when a particular upscaled model is used beyond its range of validity to provide additional insight.


Introduction
Fractures in porous media are encountered in several contexts including petroleum sector, water contamination, and nuclear waste disposal. A fracture is a thin but long domain embedded in a porous matrix and it is one of the characteristics of the heterogeneities in the medium. The flow and transport properties such as permeabilities is often drastically different from that of the matrix. The influence of these fractures on the flow behaviour is quite strong [1]. For instance, if there is a fracture network, the entire flow may take place through the fractures. At the same time, small aperture to length ratio makes it difficult to resolve the flow behaviour explicitly through brute force computations. An exhaustive survey of literature dealing with fractures is quite ambitious and for a recent survey we refer to [7].
There are several approaches to incorporate fractures in the flow models. First, to resolve them fully by considering them as equidimensional domain as the porous medium. Secondly, we may consider the fractures to be lower dimensional geometric objects embedded in the porous matrix. This will imply that the model equations are defined in heterogeneous domains with partial differential equations defined in both the porous matrix as well as on the surface of the fracture. Thirdly, we can simplify the impact of the fracture by incorporating their effect by suitably modifying the permeability of the grid cells. This paper is concerned with the second approach: to incorporate fractures as lower dimensional geometric objects embedded in the porous matrix. In practice, one combines the second and third approaches as the subsurface is full of fractures of different shapes and sizes. The smaller fractures are accounted for by suitably modifying the permeability fields whereas the larger fractures need to be treated explicitly. Thus, our approach here assumes that the smaller fractures are already incorporated by a suitably altered but given permeability field in the matrix.
The particularity in our derivation of reduced models is considering the polymer enhanced oil recovery (EOR) as a flow model. Polymer injection is a well known enhanced oil recovery technique where we mix large molecules into water to reduce its mobility and hence improve the sweep efficiency. However, in a fractured domain, the flow behaviour is strongly influenced by the fractures. Moreover, the fractures, due to their significantly different permeability and geometry, cause a large variation in the shear stresses and the velocity of the polymer [8]. Therefore, the effects of polymers in a fractured media becomes more complex. Such complexity needs to be accounted for in any polymer enhanced oil recovery model in a fractured porous medium. For example, if the flow takes place mainly along the fractured network, the sweeping efficiency will be reduced causing poor recovery. EOR applications using polymers have been widely described in literature and we refer to the works of [9,25,35,37,43] for more information.
The fracture flow models treating the fractures as surfaces have been widely used [7]. Single phase and multiphase flow models are considered in [6,23,[28][29][30]33], models with a network of intersecting fractures are in [15,17,18], transport models are discussed in [32]. Extensions to multiphysics including geomechanical effects are considered in [12,22]. Also, see references within. The advantage of treating fractures as a lower dimensional geometric object is that we no longer need to resolve the fractures through a fine meshing even while explicitly including the their effect. Resolving them explicitly as a equidimensional geometric object implies high computational costs. For example, ratio between a characteristic fracture width to a typical porous matrix is in orders of 10 −6 , resolving which using a finer mesh in the fractures will be quite expensive. Moreover, we may not be interested in the details of the pressure behaviour in the thin width of the fractures. Treating them as surfaces alleviates these difficulties and often provides us the freedom of choosing suitable and independent discretization on the surface. This approach, with some technical differences, is also known as Discrete Fracture Networks, Hybrid Fractured Models, Mixed dimensional models etc.
The reduced flow model depends on the permeability of the fracture. This is easy to motivate. Let us consider a porous domain where the matrix consists of two isolated subdomains separated by a fracture. If the fracture permeability is extremely low then it acts as a barrier. In this case, the two subdomains are decoupled. On the other hand, when the fracture has sufficiently high permeability, there would be no pressure difference in the fracture. Thus, it would be natural to define a flow model on the fracture surface that is coupled to the porous matrix. Our approach here is to identify a small parameter which is the ratio of the width of the aperture to the length of the fracture. We consider the fracture permeability as a diagonal tensor: We assume here that the fracture is filled with porous material. This justifies using Darcy scale model inside the fracture. Moreover, the filling material may also be anisotropic and for this reason we consider the permeability to be a diagonal tensor. Anisotropic permeability is well used in literature, for instance in the works of [36,38] The upscaled models are obtained by taking the limit tending to zero. We identify upscaled models for different scalings given by α and β. This should be interpreted as follows: let us say, we are given the permeability of the fracture. The question arises, which upscaled model shall we use for the fracture? The answer depends on the scale of the permeability compared to the width of the fracture. Using given width and permeability, we can compute α and β for the fracture. Once, they are known, we can determine the appropriate flow model to be used for the fracture surface. The effective models in this paper are found using a formal asymptotic approach where we expand the solution variables in terms of . We define the reduced model as the leading order term of the solution and identify the equations that are satisfied by the leading order term. Since this approach considers an expansion ansatz, we term this as a formal approach. The vanishing limit of should be considered as a mathematical exercise as in practice always remains small but positive. However, the limit equation approximates the model with the approximation getting better as becomes smaller. We also mention here that the terms reduced model, upscaled model, and effective model are synonymous and will be used interchangeably.
We make a distinction in the approaches involved in derivation of fracture models in terms of the ones sustained by rigorous mathematical analysis and the other ones by formal arguments. For a rigorous upscaling of fractured flow model, we refer to the works of [29,30,40] for a single phase flow. In case of the Richards equation we refer to [27] where mathematically rigorous convergence results are obtained for a certain range of parameters. For a study on the transport equation in similar setting we refer to the work by [32]. In [21] they use Fourier analysis to obtain coupling conditions between subdomains and obtain model error estimates when the fracture is represented as a hypersurface embedded in the surrounded rock matrix. In [19,20,31], the authors consider a thin domain with periodic coefficients for a reactive transport model and perform a rigorous two-scale homogenization to obtain interface conditions. In contrast to these rigorous works, formal arguments are used to derive the effective models in [28] which is a standard reference for such models. We also refer to the work of [3] where they study a reduced fracture model for two phase flow and the work of [18] where they study a reduced model for flow and transport in fractured porous media. A somewhat different approach, inspired from asymptotic homogenization uses asymptotic expansion of solution variables for deriving a catalogue of models for a range of coefficients. We refer to [24] and [32] for such a formal upscaling and numerical validation of the Richards equation modelling unsaturated flows and reactive transport in fractured porous media, respectively. The current work follows the approach in [32] and [24] to derive the upscaled models.
The upscaled models are typically mixed dimensional partial differential equations with the flow equation in the porous matrix coupled to flow equations on the fracture surface. Since the resulting systems are nonlinear, the discretization yields a nonlinear system of algebraic equations and we need to use nonlinear solvers such as Newton method. Since the medium consists of different homogeneous blocks, domain decomposition methods are quite appropriate. In our case, the mixed dimensional models have been solved using a domain decomposition type approach where we solve the flow in the matrix and the flow in the fracture independently and then iterate between the two to achieve the converged solutions.

Structure of paper
In Section 1.2 we describe some of the polymer models used in practice along with the Carreau model which is the one we will use later. In Section 2 we describe the problem of polymer flow in a fractured medium and in Section 3 we present the different upscaled models. In Section 4 the effective models are derived and the results from the numerical computations that illustrate the quality of our upscaled models are presented in Section 5.

Polymer model
The particularity of our derivation is considering polymer enhanced oil recovery (EOR) models. The polymer EOR technique employs injecting polymers to alter the viscosity of the injecting fluid to improve the sweep efficiency. The flow control is therefore through a suitable alteration in the viscosity. The viscosity becomes a non-linear function of the shear rate,η. In the literature, such a relationship is often expressed by the following experimentally fitted expression first given by Chauveteau (1882) where K is the permeability and φ the porosity of the porous medium. The term c is a correction factor that depends on the porous medium and the properties of the polymer, and different values may be obtained for different scalings [44]. The term u is the Darcy velocity given by Since c, K and φ are given functions of space, we can let the viscosity just be a function of velocity itself, μ(u).
This identification will simplify our presentation. The most widely used model to describe such a polymer shear viscosity relationship is the power law model, also called the de-Waele Ostwald relation [8] μ = kη n−1 (4) where k is a constant equivalent to Newtonian viscosity when the power law constant n = 1. For the case when n > 1, this law describes shear thickening behaviour, which is an increase in viscosity as the shear rate increases. Correspondingly, it describes the shear thinning behaviour, which is a decrease in viscosity as shear rate increases for the case n < 1. This power law function is strictly monotone and therefore, is unrealistic for large values of shear strain where viscosity approaches a constant value μ ∞ . This interval for the constant value μ ∞ is known as a plateau. There is a similar plateau for low shear rate where we define the viscosity as μ 0 . To incorporate the upper and lower Newtonian plateaus for μ 0 and μ ∞ an alternative to the power law model is widely known as the Carreau Model [13].
Where λ is the polymer relaxation constant,η is the shear rate and n is a power law constant (see Eq. 4). When n < 1, this model describes shear thinning with μ ∞ being the viscosity at infinite shear rate and μ 0 the viscosity at zero shear rate. When n = 1 the fluid is Newtonian, and for the case n > 1 the model behaves as a power law fluid describing shear thickening. One of the drawbacks of the Carreau model is that it cannot describe both shear thinning and shear thickening for the same fluid [4]. In practice, complex compositions of fluid mixtures are used that show both shear thinning and shear thickening behaviour in the same fluid at different locations or time. There have been suggestions in literature, see e.g., [14] to include both shear thickening and thinning behaviour in the same fluid by including more parameters. One of the models proposed is where τ r is the polymer relaxation time, λ 2 is a second time constant and n 2 a second power law exponent. This extended Carreau model is also called the Universal Viscosity Model (UVM) and it describes the apparent viscosity as a function of 10 parameters.
In what follows, we study the Carreau model and present the derivation of reduced models (Fig. 1). with the fracture domain f as a rectangle symmetrically embedded. The domain and the subdomains are given by: and for the limit → 0 we denote Here, the domain contains a single fracture f which divides into two subdomains 1 and 2 . We denote the boundary of 1 2 , and f by 1 2 , and f , respectively. The interfaces between the subdomains are denoted by 1 and 2 . Here 0 denotes the width of the fracture and is much smaller than the length and hence, is a small parameter. Extending this to an open, bounded, convex polygonal domain containing a non-self-intersecting fracture is straightforward.
As already indicated above, resolving the thin domain f requires a meshing that is of the order of which is typically very small. This will be computationally costly. The idea is therefore to perform the dimensional reduction of the fracture model and obtain an equation on the interface such that the solution of the dimensionally reduced model approximates that of the original problem. Naturally, the processes in the fractures are coupled to the porous subdomains and hence the resulting reduced models will be of a hybrid-dimensional (or mixed-dimensional). As we will see depending on the fracture hydraulic properties, there is a catalogue of models that we obtain.
We consider a non-Newtonian single phase flow in governed by Darcy's law and the mass conservation equation.
Here, μ(u ) is the viscosity as a function of velocity, and it can be described by the extended Carreau model in Eq. 6 together with the relationship between shear stress and velocity described in Eq. 2. q is a source/sink term, and the permeability in j , j = 1, 2, and f is given in Eq. 10 For simplicity we assume that the permeabilities K 1 and K 2 in 1 and 2 equal 1 in both x and y-direction. Even though the permeabilities in 1 and 2 are constant, we denote it with the -superscript for notiational convenience. In the fracture domain the ratio between the length and the width motivates us to describe the permeability as an exponent of fracture width Here, α, β are given real numbers and are the parameters that describe the permeability in the fracture. We will see that as we vary the parameters α, β, different reduced models are obtained. In practice, they are estimated as follows. We assume that the fracture is given along with its permeability and geometric dimensions. The latter provides us which is fixed for this given fracture. Next, we compute the scaling of the permeability in both x and y directions with respect to . This yields the parameters α and β.
For the sake of convenience homogeneous Dirichlet boundary conditions are considered. The interface conditions are the continuity of the pressure and the continuity of normal flux. We get the complete set of equations for our fracture problem where j = 1, 2 and i = 1, 2, f and K f and K j is defined in Eq. 10. In the above system of equations, we have introduced as a superscript to emphasize the dependency of the solution variables on . In these models it is assumed that the pressure is continuous at the interfaces separating the matrix blocks and the fracture. We mention other works such as [41] where an extended pressure condition is used instead of the continuity of pressures at the interface. Such conditions can explain the trapping of oil when the flow moves from a high permeable region to a lower one [42]. In case of oil recovery, such effects may be important. Homogenization approaches that derive such models are in [34,39]. However, such models are posed directly at a finite model (where the fracture and matrix are equidimensional). In our case, we have considered a simplified approach by assuming pressure continuity. However, these extensions can be obtained in the framework that we present by suitably altering the starting microscopic (fixed ) model. This should be contrasted with the upscaled models obtained here where we do get discontinuity in the pressure across the interface as well. However, this is in the effective equations rather than in the model itself, see also [2,5,11,24]. In Section 3, we take the above model, Eq. 11, as the starting one and study the limit of tending to zero. Accordingly, we perform a formal asymptotics, where we expand the solution variables in powers of . We identify the equations satisfied by the leading terms and define them to be the reduced models. Naturally, they depend on the hydraulic properties of the fracture and in our description on α, β. Therefore, we get a catalogue of reduced models depending on the values of α and β.

Catalogue of reduced models
We state the main results in this section. We begin by classifying our effective models in three types. These are: i) models with continuous pressure across the fracture-matrix interface, ii) models with discontinuous pressure across the interface and iii) a model with disconnected pressure across the fracture-matrix interface. The difference between the last two is the pressure distribution in the effective model inside the fracture. In case ii), the pressure is given by a differential equation inside the fracture which is coupled to those in the matrix while in case iii), the fracture acts as no-flow boundaries for the flow model in the matrix (Sections 3.2 and 3.3). The derivation of the reduced models is done in Section 4. For all the models [u j ] denotes the jump in the flux [u j ] = u 1 · n − u 2 · n with n being the outward normal on pointing towards 2 .

β < 1: Model with continuous pressure in the fracture
When β < 1; we get models that have pressure continuity across the fracture-matrix interface, independent of α. The permeability in y-direction in the fracture allows the equilibration of the pressure in the y-direction so that the fracture collapses on a surface and takes the value of the trace of the pressure from the matrix side.

Effective model 1
For α > β − 2, β < 1, j = 1, 2 we get the following model In this case, the value of permeability in the y-direction is sufficiently high so that the pressure is equilibrated in this direction. Accordingly, the pressure becomes uniform along this direction. Moreover, this case also implies a low permeability in the x− direction so that the resulting model has no flow along this direction. This makes the effective problem as if the fracture is an interface playing no role, and the pressure and the normal fluxes in the matrices are continuous. Accordingly, the flow equation in the matrix can be solved independently without considering the flow equation in the fracture. The two interface conditions for the matrix flow equations are the continuity of fluxes and that of the pressure. Once the flow equation is solved for the matrix, fracture pressure can be computed by using the continuity of pressure to the trace of the matrix pressure on fracture surface. Similarly, the flux can be computed by using the pressure in the fracture.

Effective model 2
Effective model 2 As in the previous case, we see that the permeability in the y− direction is large enough for the equilibration of the pressure there. Moreover, the elliptic equation in the fracture retains its character though now defined on the collapsed interface. The continuity in the pressure is retained through the traces of the pressure from the matrix side. However, there is still a jump in the flux that may be non-zero which is on the right hand side in the fracture equation.

Effective model 3
For α < β − 2, β < 1, j = 1, 2, the effective model takes the form The pressure as before is independent of the y variable. Moreover, the pressure continuity across the fracture matrix interface is retained in the effective model. Furthermore, the pressure in the x-direction also gets equilibrated as the permeability in the x-direction is quite large. This means that the pressure in the fracture is determined by solving an ordinary differential equation defined on the fracture surface. This pressure acts as the boundary condition for the matrix equation with homogeneous boundary conditions, p f = 0. In case of homogeneous Neumann boundary condition on the boundary of the fracture, the pressure on the fracture surface becomes constant. This constant is determined by the condition [u j ] = 0. As the numerical tests presented later show, the mass conservation is retained in this case.

β = 1: Model with discontinuous pressure
If the y-direction permeability is proportional to the fracture width, that is, β = 1, effective model results in which the pressure is discontinuous across the interface fracturematrix interface and there is a jump in pressure over the fracture. In this case, we realize that the fracture cannot be collapsed on an interface. Thus, the usual sense in which we obtain an effective equation defined on an interface, that is by averaging across the fracture aperture cannot be justified. Hence, we get a two scale model where the pressure jump is determined by an elliptic differential equation over the fracture domain f := (− 2 , 2 ) × (0, 1). In this case, it is natural to rescale the fracture domain so that it is independent of the value of . We define new variables in the rescaled domain. ξ = x and η = y . The fracture domain in the new coordinate system becomes f := (− 1 2 , 1 2 ) × (0, L). For notational convenience, we define,

Effective model 4
For α > −1, β = 1, j = 1, 2 gives the following model Effective model 4 As for the other effective models, the flow equations in the matrix remain unchanged. Here, e 2 is the unitvector in η-direction. We see that the dominant flow is going in the η direction. For any fixed ξ , p f is a solution of a second order boundary value problem in the vertical direction with the unknown p f being a function of η with the boundary conditions at η = − 1 2 and η = 1 2 . These boundary conditions are given by the traces of the pressure from the matrix. Since the traces of the matrix pressure are themselves unknowns, the flux conditions at the boundaries η = − 1 2 and η = 1 2 complete the system. There is no horizontal coupling as the ξ variable enters the model as a parameter. In other words, the flux has only one component. In the case here we have that the flux is constant, however, in cases with external sources this would not be true.

Effective model 5
Effective model 5 This scaling of the permeability in the fracture domain is the critical one. In this case, the permeability in the ξ and η directions after the fracture domain is rescaled is such that it gives the same flow equation in the fracture as in the matrix. This may include a jump in the flux as well as the jump in the pressure across the fracture surface if one views it with respect to the matrix flow equations. These are determined by a differential equation in the rescaled fracture domain f . We see the emergence of a two-scale nature of the upscaled equation.

Effective model 6
For α < −1, β = 1, the model takes the form Effective model 6 Here, the diffusion in the ξ direction is much stronger so that the other terms in the fracture flow model vanish in the limit. This makes the pressure p f satisfying a boundary value problem involving an ordinary differential equation. This explains the derivatives in f being an ordinary derivative in contrast to the partial derivatives in porous matrix.
Together with the boundary conditions on ξ = 0 or ξ = L, the pressure in the fracture can be obtained. This acts as Dirichlet boundary condition for the matrix flow equation. Note that in the simple case of taking homogeneous boundary conditions at ξ = 0 or ξ = L, the pressure in the fracture would be zero. In case of homogeneous Neumann boundary condition on the boundary of the fracture, the pressure in the fracture subdomain becomes constant. This constant is determined by the condition [u j ] = 0.

β > 1: Model with disconnected pressure
For the case of β > 1, the permeability in the fracture is much smaller than that of matrix. In this case, the fracture becomes a non-conducting barrier and the flow is confined to the matrix. Accordingly, we get a decoupled model where the two solid matrix domains 1 and 2 are entirely decoupled from each other and can be solved separately.

Effective model 7
We have the following model for j = 1, 2.
div u j = q j in j The last two equations are the boundary condition for the matrix flow and shows that the fracture interface acts as a no flux condition. The "fracture" thus acts as an impervious barrier. Under these conditions, the fracture surface provides sufficient number of boundary condition for the flow on the either side of the matrix. The matrix flow equations can be solved uniquely. Accordingly, the equation in the fracture accordingly becomes irrelevant.

Remarks on the upscaled models
Observe that unlike other effective models, in the Effective models 4, 5 and 6 the solution in the fracture depends on η as well. Also, η = − 1 2 in the fracture corresponds to the top boundary of the matrix domain 2 (at y = 0−). Similarly, 1 has the bottom interface at y = 0+, connected to the fracture domain boundary at η = 1 2 . This is a "two-scale" type of model where the permeability has the scaling such that the dimensionality reduction of the fracture interface is not justified. Instead, the original problem posed in a larger domain is reduced to three coupled sub-problems. Similar models, in the context of reactive flow, are derived in [31]. This is in contrast to the Effective models 1-3 where the solution in the fracture domain is independent of η and the fracture collapses as an interface. Also the two interfaces from the porous matrix sides coincide. We make a brief comparison with the fracture models that are widely used in practice (see e.g., [7]). We refer to [28] where similar models are derived for single phase flow, to [3] for two phase flow models and to the recent work of [21] where Fourier analysis is used to derive coupling conditions between subdomains. The key difference is in the methodology of derivation and the starting assumptions. Whereas in the references cited here, the authors make an assumption on the pressure profile inside the fracture, our starting assumption is the scaling of permeability of the fracture based on . However, once the assumption on scaling of the coefficients is made, the derivation follows without requiring any assumption on the solution. In contrast, the derivation in [28] and in [3] takes place by first integrating the flow equation along the transverse direction in the fracture subdomain. Using the continuity of fluxes at the matrix/fracture interfaces, this yields a surface equation with the jump in the matrix flux term as a source term. Further, a closure relationship is postulated for the pressure profile in the fracture. In the work of [21] the derivation of coupling conditions relies on a Fourier transform of the physical unknowns in direction tangential to the fracture and an elimination of the fracture unknowns in Fourier space. In contrast, our approach does not postulate any closure condition on the pressure inside the fracture as this is part of the solution. Also we mention that in the references cited, the closure condition introduces a parameter in the effective model for the fracture. Here, we have a catalogue of models and no additional parameter is necessary. Moreover, for some of the regimes considered here, for nonlinear models such as the Richards equation, the derivation is sustained by mathematically rigorous proofs (see [27]). We also refer to [21] for rigorous modelling error estimates for models in the spirit of [28]. For more detailed comparison between the existing models in literature and the current approach, we refer to [24].
In this paper a simple geometry for the fracture embedded in a porous matrix and a simplified polymer flow model having natural interface conditions across the matrix/fracture interface are considered. These interface conditions are the continuity of pressure and those of the normal fluxes. The formal upscaling procedure that we employ here allows us to consider more complicated flow models involving more physics or complex constitutive models. Our objective here is to illustrate the role of properties of the fracture (compared to the ones in the adjacent blocks) in determining the appropriate Discrete Fracture Network (DFN) model for polymer EOR models. For further uses, we mention that this approach can be used to provide arguments for the mixed dimensional models proposed in [10] based on physical arguments. In the same spirit, the reduced models are derived for fracture networks in [16] building on the work of [28]. Our work can be used to similarly provide a hierarchy of reduced models based on the properties of the fractured networks. Our formal approach can be adapted to several extensions. Matrix heterogeneity can be immediately included as it does not change any of the descriptions below. Hence, the findings here can be incorporated to obtain other types of flows such as compressible, multiphase flows, or reactive flows or other broader models.

Derivation of upscaled models
We start with the subdomain equations followed by the fracture equation. Together with the interface conditions, we will be in a position to derive the upscaled models for different values of α and β.

The subdomains 1 and 2
The derivation of the upscaled model in the subdomains 1 and 2 is rather straightforward. Our starting model is We first note that as vanishes, the domains j tends to j . Moreover, as stated above, we have chosen the permeability K j = 1. This is not restrictive as extension to any positive definite anisotropic and diagonal tensor is straightforward. The viscosity μ is a function of velocity u j as already explained in Eq. 5. The dependence of the viscosity on velocity u j makes this equation non-linear.
For simplicity, we now assume that the source term q j = q j is independent of and of order O(1). In case of a two-dimensional geometry, we can explicitly write the derivatives in the Eq. 12 for the respective domains j . Now we make a homogenization ansatz. We assume that the pressure can be expanded as and the same can be done for velocity To take care of the viscosity, we assume it to be a smooth function. The Taylor expansion then provides the expansion: where μ (u j, 0 ) is the derivative of μ when u = u j,0 . We can find the effective equation by neglecting the higher order terms and insert the lowest order term into Eq. 13. The effective equation in 1 becomes Accordingly, the effective equation in 2

The fracture f
The fracture flow model gives, Here For simplicity, we set q f = 0. We eliminate the flux above to have an elliptic equation for pressure, The fracture domain is dependent on and a simple rescaling maps it to an − independent domain. Indeed, define ξ = x, η = y to rewrite the equation for the fracture flow as In the above equations, we have retained the notation to denote the same symbol for pressure as in the original domain, that is, p f (ξ, η) = p f (x(ξ ), y(η)). As in the case of subdomains, we have the homogenization ansatz: For the viscosity, again we use the Taylor expansion where μ denotes the derivative of the viscosity μ at u f,0 .
The equation for f can be written as

Interface conditions
Along the interfaces 1 = 1 ∩ f and 2 = 2 ∩ f we assume that there is a continuity of normal fluxes and pressures across the fracture/matrix interfaces. For j = 1, 2 Using the expansion for the pressures, we obtain at In terms of η = y , we have at j , j = 1, 2 Equating the respective powers of for the first equation, and note that as → 0, j converges to j , we have for j , j = 1, 2 P f,0 = P j,0 , P f,1 = P j,1 .
For the second equation in Eq. 26, we consider different values of β. For β > 1, we have upto ), For β = 1, For β < 1, we keep it in the same form and use it in later derivation.
It is clear from the preceding discussions that the interface conditions at the fracture/matrix interface depends on the permeability scaling; accordingly, the effective equations will vary depending on the permeabilities in the x− and y− directions. We study this for each of the cases below.

Case 1 β < 1
Let us first note an energy inequality that will be used in this part. The following estimate is easily obtained using standard variational techniques. The above inequality is obtained by testing the flow (11) with p f , p 1 , p 2 and integrating over the respective domains f 1 , and 2 . Performing partial integration and using interface conditions provide us the above estimate. Here, C is a generic constant that depends on the boundary condition and the right hand side. Using the bound that the viscosity is bounded by below μ f ≥ μ 0 , we get Switching to the ξ, η co-ordinates with ξ = x, η = y transforms the fracture domain f to f , we get Here, we have reused the same notation for p f but now as a function of ξ, η. This allows us to conclude for β < 1, ∂ η p f = 1−β ) and the leading order terms, P f,0 is independent of η.
In this case, the first term is dominated by the second. Hence, Eq. (21) is reduced to However, as we already concluded P f,0 is independent of η.
In this case, the fracture therefore simply disappears. That is, across this interface, the flux continuity and the traces of pressure from the subdomains become equal to that of the fracture. There is no flow along the fracture.
Using the interface condition, we get d dξ Here, [− 1 μ(u 0 ) ∇P ] = − 1 μ(u 0 ) ∇P 1,0 · n + + 1 μ(u 0 ) ∇P 2,0 · n − with n + and n − denoting the outward normal on 1 and 2 respectively. Using P f,0 is independent of η, and reverting back to x, y coordinate, we obtain Together with the equality of traces of the pressures P f,0 = P 1,0 = P 2,0 on . This completes the model. Note here that the fracture model has a partial differential equation defined on its collapsed surface with the jump of the fluxes from the matrix in the right hand side. This is the most commonly model used in the literature. This corresponds to the situation when there is a flow taking place along the fracture surface with the flow from the matrix side entering or exiting. In this case, the flow equations are fully coupled: the matrix flow and the fracture flow equations together with the equality of traces of the pressures determine the pressure profile.

Subcase 1.3 α < β − 2 In this case, the first term dominates and the equation takes the form
Here recall is the collapsed surface of the fracture. The pressures on the fracture surface can be found by solving this ordinary differential equations and using the boundary conditions at the ends of the . In the case of homogeneous Dirichlet boundary conditions, this would lead to zero solution. This can be then used as boundary condition for the matrix flow. As we show in the numerical examples, the model works well also for non-homogeneous Dirichlet boundary conditions.

Case 2. β = 1
We again begin from Eq. 21 We subdivide this in three regimes.

Subcase 2.1
For α > −1, β = 1 we multiply the above equation, Eq. 21, by to make the second term in the left hand side free of and noting that the first term on the left has α+1 as the coefficient which vanishes as tends to zero. Thus, we have for the leading order term, The above holds for every ξ . The boundary conditions for this ordinary differential equation are obtained by the interface conditions at η = ± 1 2 . Thus, we have, for each ξ (and recall x = ξ ), P f,0 (ξ, η = −1/2) = P 2,0 (x, y = 0).
Note that P 1,0 and P 2,0 are unknowns satisfying the polymer flow equations in the subdomains 1 and 2 . The above ordinary differential equation thus cannot be solved unless we specify further conditions. Indeed, the continuity of normal fluxes is retained via the interface conditions ∂ η P f,0 ξ, η = − 1 2 = ∂ y P 2,0 (x, y = 0). (37) Note that P f,0 = P f,0 (ξ, η) and η = ± 1 2 corresponds to the top/bottom boundary of f . At the same time, this coincides with the boundary y = 0 of the subdomains 1 and 2 . Moreover, the equations in all the subdomains are coupled together. This case therefore brings out the "twoscale" nature of the effective equations as alluded before. Here, the equation for the fracture is set up in a domain f and cannot be collapsed on an interface as are the cases for β < 1 as shown in the next section. Further, the pressure at the matrix fracture interface as seen from the subdomains 1 2 are discontinuous. At the same time, we also have a well-defined pressure equation in the fracture that couples the traces of the two subdomain pressures and allows these traces to be discontinuous. This is in contrast to the case β > 1 where the traces of the subdomain pressures may be discontinuous as well, however, the fracture equation becomes irrelevant as it behaves as an impermeable barrier. The present case also shows that in the regime discussed here, the averaged model needs to consider the full details in the fracture flow.

Subcase 2.2
For α = −1, β = 1, the exponents of balances for both the terms on the left hand side. This implies after multiplying the equation on both sides by , we obtain The solution requires interface conditions. These conditions are the continuity of the pressures and that of the normal fluxes The continuity of normal fluxes is given by ∂ η P f,0 (ξ, η = − 1 2 ) = ∂ y P 2,0 (x, y = 0).
The same comment as above for this model being twoscaled holds here as well. Moreover, as in the previous case, the effective equation cannot be collapsed on an interface. In the list of effective models, together with the subdomain equations, this is Effective Model 5.

Subcase 2.3
For α < −1, β = 1, we follow an analogous argument and obtain d dξ We get an ordinary differential equation which can be solved and is independent of η (independent of y-direction).
Together with the subdomain equations, this gives the Effective Model 6.

Case 3. β > 1
In this case, the interface conditions simplify and yield equations such that the two subdomains 1 and 2 becomes decoupled and obtain no-flow interface conditions. This is as vanishes, for β > 1, β−1 → 0 and since other terms are assumed to be of O(1), we obtain 1 μ(u 0 ) ∂ y P 2,0 = β μ(u 0 ) ∂ y P f,0 = 0 at 2 (45) As the interface condition for the subdomain takes the form of a no-flow (homogeneous Neuman) boundary condition, the subdomain equations can be solved independently. We have thus, a simplified case when the Darcy equations in the subdomain are retained with the fracture/matrix interface acting as a no-flow boundary condition.
We call this the interface conditions from the Effective Model 7. Since the model is complete with the no-flow boundary condition for the subdomains, the pressure inside the fracture becomes irrelevant. Thus, irrespective of the value of α. the boundary conditions gives us that there is no flow in y-direction at the boundaries 1 and 2 . This case represents the impermeable barrier in the effective model. We will validate this model in the numerical computations later.

Summary in a table
We summarize the above in the table below.

Numerical results
In this section we validate the results of the upscaled models numerically. As the derivation shows, the upscaled models are valid for a certain range of α and β. Thus, the validation exercise here is to compare the effectiveness of a particular upscaled model in its range of validity and the errors caused when it is used outside this range.
Our reference geometry as defined in Eq. 7 and for the reduced model the geometry is as defined in Eq. 8. For the discretization, subdomains 1 and 2 are both split up into 100 × 50 grid cells. The edges of the domain have been given Dirichlet boundary conditions. Since we have a variety of upscaled models, solving them requires different approaches. The simplest is the case of no-flow boundary condition (effective model 7) where we prescribe a no-flow boundary condition for the interface between the two subdomains. In the case of "two-scale" model (effective model 4, 5, and 6) we define an iterative procedure. While solving the effective model in subdomain f , we use the pressure boundary condition from the two subdomains. Whereas, while computing the solutions in subdomains 1 and 2 , we take the flux computed from the subdomain f and use it as a boundary condition. In the effective models 1, 2, and 3, again we perform an iterative procedure. We solve the subdomain problems 1 and 2 using the pressure boundary conditions computed by solving the interface problem. Moreover, the flux is computed at the subdomain/fracture interface and the jump in the fluxes is added to the right hand side in the fracture interface model. Our implementation has been done in the Matlab based open source code MRST [26].
For all the simulations we have used a finite volume scheme on a static uniform grid with rectangular cells. The fluxes are implemented first with a two point flux approximation and the pressures are defined at the center of the cells. Note that the Carreau model is a nonlinear flow model, which is solved by an iterative procedure. The iterative scheme consists of using the viscosity computed from the iteration-lagged flux in the flow problem, hence linearizing the equations. After solving them, we compute the flux and calculate the updated shear stress, which is used in the Carreau model equation to provide the updated viscosity. This process is repeated until convergence. The test for convergence is the difference of pressure between two successive iterations. The norm used is L 2 and the tolerance is taken as 10 −6 . For all our cases we have used the standard Carreau model, and the parameters used are μ 0 = 10 cP (1 × 10 −2 Pa.s), μ ∞ = 1 cP (1 × 10 −3 Pa.s), c = 1, λ = 1 s, and φ = 1.
Our approach for the validation is as follows. We select certain values of which determines the reference geometry and the permeabilities. We compute the pressure for the full -model resolving the fracture thickness. i.e. we solve the pressure equations (13) - (19) with a finite value of . This forms our reference solution. Corresponding to the chosen value of α and β, we have the recommended upscaled model and we compare it to the full -model. An example is shown in Fig. 2. Here Fig. 2a gives the pressure in the domain for effective model 1, and Fig. 2b gives the pressure when the full -model is used. In both the Fig. 2a and b, we have fixed α = 2 and β = −2 and = 3 50 . In the next sections we study in further details the difference of solutions for the effective models and themodel. We plot the pressures at the cells following three lines parallel to the y-axis at x = 0.4, x = 1 and x = 1.6, and also plot the difference between the two pressures.

The average model: effective models 1-3
We investigate when the fracture has sufficient permeability in y-direction so that the pressure is constant in that direction. We choose α = 2 and β = −2 which falls in this regime. The appropriate effective model for this case is Effective Model 1. We solve both the reference and the effective model in this case. For both the problems, we discretize the domain having 100 × 100 grid cells, and we choose the fracture width to be = 3/50. The boundary conditions are p = 0 at top (y = 2) and p = 1 at bottom (y = 0), and at the sides we have p = 0 at x = 0 and p = 1 at y = 1. Note that the boundary conditions may affect how well the effective model approximates the pressure. This boundary condition is different from the model boundary condition (homogeneous Dirichlet boundary condition) that we chose for the derivation, however, the boundary condition chosen here is more interesting for pressure profile. The pressure for the Newtonian flow is plotted in Fig. 2 and again in Fig. 3a for fixed values of x at x = 0.4, x = 1 and x = 1.6. The difference in pressure between the average model and the -model is plotted in Fig. 3b. As we see the two models are overlapping and the difference in pressure is of the order of 1e − 7. We conclude that the average model is a good effective model in this case. Now we take the same simulations for the non-Newtonian case which is of interest in this paper. In Fig. 4 we have used the average model for the non-Newtonian flow. The only thing that has changed from the plot in Fig. 3 is that the n parameter in the Carreau is now equal to 0.5, and this makes the fluid non-Newtonian. As we see from the Figure, the average model is still a good effective equation.
We see above that the average model works fine for the case where α = 2 and β = −2. This shows that the effective model works fine in its range of validity. To show how different boundary conditions affect the results we have plotted the pressure in Fig. 5a for a case where the boundary conditions are p = 1 at y = 0 and p = 0 on all other sides. The difference between our upscaled model and reference model can be seen in Fig. 5b. To show that the upscaled model is even better for small grid sizes we have plotted the pressure for a case where the domain is divided into 1000 × 1000 grid cells and where = 3/1000. The approximation becomes better with more grid cells. Note also that the error here is of the order of 1e − 7 which is of the same order as the tolerance of the non-linear solver (Fig. 6). Now, we test the average model for a wide range of different α and β values including outside the range of validity of the model. This will give us the performance of any particular model and its use without regard for the permeability values of the fracture. In Fig. 7 we have found the relative difference between the average model and the reference model for different α and β and plotted the difference in the figures. In Fig. 7a we have the range of validity for a Newtonian flow, and in Fig. 7b we see the validity range for the average model for a non-Newtonian flow. As expected in both cases, Newtonian and non-Newtonian, the average model works fine for β < 1. At the same time, the approximation becomes poor if the model is  pushed to be used for permeability ranges outside its validity range. This shows that a single upscaled model cannot be used for arbitrary values of the fracture permeability.
We make a remark on the reduced model 3. Here, the pressure in the fracture is determined by an ODE defined on the fracture interface and the pressure obtained thereby is used as Dirichlet boundary condition for the matrix side. This raises the issues of conservation of mass in the reduced model. In this case, the permeability in the x− direction is quite strong which dominates the flow behaviour. We distinguish the two cases for the boundary conditions for the fracture subdomain. In the homogenenous Neumann boundary condition, the pressure becomes constant on the fracture surface. However, this constant is determined by employing the mass continuity condition   We split our domain into 1000 × 1000 grid cells. = 3/1000. This is the case where we have a non-Newtonian case. The n in the Carreau model is equal to n = 0.5 Fig. 7 The domain error gives us that the validity of the different models are dependent on both the scalings α and β  equal to 4 and is then used as Dirichlet boundary condition for the matrix side. We plot the flux jump in both the cases in Fig. 8. Moreover, for different values of , we compute the error in the mass balance in the two cases in the maximum norm. We tabulate the results in Table 1.
Furthermore, in the homogeneous Neumann case, we computed the pressure profile in the reference solution for different values of pressure boundary conditions on the top and at the bottom. We verified numerically that the pressure profile is nearly constant in the fracture and assumes the value that satisfies the condition [u j ] = 0. Test case 2 for effective model 3 In this case, the boundary conditions for the reference domain are: At x = −1, p = 4; x = 1, p = 8; y = −1 − = 1; y = 1 + = 20. For the reduced model, the boundary conditions are: At x = −1, p = 4; x = 1, p = 8; y = −1, p = 1; y = 1, p = 20. We choose = 0.1, 0.05, 0.025 and compute the pressure profiles for both the reference and the reduced model. In the reduced model, the pressure profile on the fracture is given by 2x + 6 that is consistent with the two boundary conditions at x = −1, and x = 1. We then compute the jump in the flux across the fracture in the reduced model and for the reference model we compute the net flux in the fracture subdomain from the matrix side. We compare the two values in Fig. 9.
Moreover, as before, for different values of , we compute the error in the mass balance in the two cases in the maximum norm. We tabulate the results in Table 2 This suggests the convergence to the reduced model and also suggests a rate of first order convergence. These numerical tests indicate that the reduced models do not introduce any additional errors in mass balance.

Two-scale model, effective models 4-6
This covers the effective models 4, 5, and 6 for the case when β = 1. Our grid for the upscaled model becomes 100 × 150 where 100 × 50 of the cells represent the fracture, and the boundary conditions are p = 0 to the left, p = 1 to the right, p = 0 at the top and p = 1 at the bottom. The plot in Fig. 10a is again the pressure for given lines at fixed x values, and Fig. 10b is the difference between the reference model and the two scale model along the same lines. From the plot it is obvious that the two scale model is a good effective model for the case when α = 1 and β = 1 for the Newtonian flow as expected from above calculations.
In Fig. 11 we again show the same results for a non-Newtonian flow, and conclude that the model is a good effective model for the case when α = 1 and β = 1. It is worth to notice that the polymer have led to a higher jump in flux over the fracture for this case.
The two-scale model has been used to solve for the pressure for a range of α and β values, and the relative error for each case is plotted in Fig. 12. Here we see the Newtonian case in Fig. 12a and the non-Newtonian case in Fig. 12b. Both the plot for Newtonian flow and non-Newtonian flow show that the two-scale model is a good effective model for the case where β = 1. From the two plots it is clear that the polymer behaviour also impacts the accuracy of the approximation. Moreover, as already stated above in the case of average model, the two-scale model    is a good approximation in its range of validity, and the approximation becomes poorer as one goes further from its validity range.

The decoupled model, effective model 7
This corresponds to the case when β > 1. In this case the fracture is much less permeable that the surrounding matrix blocks leading to a decoupled model in the effective equations. For our decoupled model we solve the two subdomains 1 and 2 separately with a no flow boundary condition at the fracture interface. The two domains 1 and 2 are made by 100 × 50 grid systems and the boundary conditions around is p = 0 to the left, p = 1 to the right, p = 0 at the top and p = 1 at the bottom. For our pressure plots and relative differnce plots in Figs. 13 and 14 we have chosen α = 2 and β = 2. From our Newtonian case we see that the difference between the effective model (decoupled model) and the reference model is of the order of 1e − 3%, and hence a good fit. For our non-Newtonian flow the behaviour is similar as shown in Fig. 14a.
The decoupled model is used at a variety of α and β values and compared to the reference model for each case. The relative error in every case is plotted in Fig. 15. Both for the Newtownian flow in Fig. 15a and non-Newtonian flow in Fig. 15b, and as expected from the previous calculations the decoupled model is a good effective model for the cases where β > 1. We also see from our plot that the decoupled model can be of use in some cases where β < 1 where the value of α is sufficiently high.
Finally, we comment on the solution procedure. As stated above, we have used an iterative procedure to solve the  15 The domain error gives us that the validity of the different models are dependent on both the scalings α and β flow model. This iterative process has two loops, an inner loop where the viscosity is updated for given boundary conditions, and one outer loop where the interface conditions are updated between the fracture and the surrounding domains(s iterations). To show how our models converge we have used the L 2 norm to find out how much the pressure in our domain change between each iteration of our outer loop. We calculate the average pressure difference between the domain after s iterations and compare it with the average pressure in the previous iteration. The average pressure difference is plotted in Fig. 16 for the 30 first iterations. Note that we do not plot convergence for the decoupled model as the boundary conditions do not change.

Conclusion
We consider a two-dimensional fractured porous medium where the fracture has a thickness . For the flow model, we have considered a single phase polymer enhanced oil recovery model. This model consists of a non-linear Darcy equation where the viscosity depends on the flow through the well-known Carreau type model. For this model, we perform an upscaling of the flow model which is obtained as the limit of tends to zero. The particularity in our approach is the anisotropic permeability in the fracture that is taken as the exponent of characterized by two real numbers α and β combined with polymer enhanced Fig. 16 Pressure change between iterations using the L 2 norm. For the average model we have used α = 2 and β = −2, and for the two scaled model we have used α = β = 1. Note that the y-axis is logarithmic oil recovery (EOR). We obtain a catalogue of models for different values of α and β. Using numerical computations, we study the suitability of our upscaled models in their regimes of validity. Moreover, we study the sensitivity of the upscaled models when a particular upscaled model is used beyond its range of validity. We show that any particular upscaled model, if used beyond its validity regime, performs quite poorly. The approach is a formal upscaling one that can be adapted to other models.
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/.