Modelling the Impact of Anisotropy on Hydrocarbon Production in Heterogeneous Reservoirs

Effective and optimal hydrocarbon production from heterogeneous and anisotropic reservoirs is a developing challenge in the hydrocarbon industry. While experience leads us to intuitive decisions for the production of these heterogeneous and anisotropic reservoirs, there is a lack of information concerning how hydrocarbon and water production rate and cumulative production as well as water cut and water breakthrough time depend on quantitative measures of heterogeneity and anisotropy. In this work, we have used Generic Advanced Fractal Reservoir Models (GAFRMs) to model reservoirs with controlled heterogeneity and vertical and/or horizontal anisotropy, following the approach of Al-Zainaldin et al. (Transp Porous Media 116(1):181–212, 2017). This Generic approach uses fractal mathematics which captures the spatial variability of real reservoirs at all scales. The results clearly show that some anisotropy in hydrocarbon production and water cut can occur in an isotropic heterogeneous reservoir and is caused by the chance placing of wells in high-quality reservoir rock or vice versa. However, when horizontal anisotropy is introduced into the porosity, cementation exponent and grain size (and hence also into the permeability, capillary pressure, water saturation) in the reservoir model, all measures of early stage and middle stage hydrocarbon and water production become anisotropic, with isotropic flow returning towards the end of the reservoir’s lifetime. Specifically, hydrocarbon production rate and cumulative production are increased in the direction of anisotropy, as is water cut, while the time to water breakthrough is reduced. We found no such relationship when varying vertical anisotropy because we were using vertical wells but expect there to be an effect if horizontal wells were used.


Introduction
Newly discovered hydrocarbon reservoirs are increasingly complex, heterogeneous and anisotropic (Glover et al. 2018;Caruana and Dawe 1996). Such reservoirs contain heterogeneities at different scales, exhibited over 9 orders of magnitude, from micrometres up to several kilometres (e.g. Zhou et al. 2003). Such heterogeneities can arise, for example, from local changes in deposition (Jolley and Fisher 2010), post-depositional history (Speight 2014), diagenesis (Rashid et al. 2017(Rashid et al. , 2015a and fracturing (Glover et al. 1997). In cross-bedded anisotropic reservoirs, the porosity, capillary pressure, wettability and permeability may vary significantly spatially, consequently affecting the initial distribution of fluids within the reservoir rock and their flow during production (Deng et al. 2014). Steadily improving knowledge of the processes occurring in the subsurface, together with better remote monitoring measurements and the advent of semi-automated approaches to reservoir characterisation and modelling, all require that newly developed reservoir models optimise the use of all available subsurface data, especially for complex unconventional heterogeneous and anisotropic reservoirs (Branets et al. 2009).
Since anisotropic hydrocarbon reservoirs were recognised in the literature (e.g. Young et al. 1964), there have been many studies on the effect of anisotropy on almost every rock physical property, from petrofacies (Marco et al. 2007;Johansen et al. 2004;Cosen et al. 1994), porosity (Huang et al. 2015;Zhao et al. 2013), elastic properties (Chemali et al. 1987;Nathan et al. 1990;Johnston and Christensen 1995) and permeability (Huang et al. 2017;Kwon et al. 2004, Cosan et al. 1994. In each case, the study is confined to the micro-scale processes which lead to anisotropy, rather than investigating how these properties conspire to affect the macroscopic properties of the reservoir such as changes in the local fluid saturation and pressure drawdown during production.
In the case of hydrocarbon reservoirs, the spatial distribution of rock properties at all scales, from the microscopic to basin scale, controls flow. It is the understanding of the way these properties change and how they influence reservoir-scale production that represents the major challenge and which is addressed by this work.
The first part of this challenge is the measurement of the reservoir's properties at all scales, or at least the development of an approach which can effectively extrapolate and rescale those properties. Fractal mathematics (Saupe 1988;Mandelbrot and Van Ness 1968) has a role to play in this. Since fractals are self-similar at whatever scale they are viewed, fractal ideas can be used to extend the effective range of our application of a value made at one scale to a different, higher or lower, scale. Hence, the distribution of porosity, measured at say the well-logging scale, can be fractally extrapolated up to the reservoir scale and down to a millimetric scale. This is the idea that drove Al-Zainaldin et al. (2017) to develop what was called Advanced Fractal Reservoir Models (AFRMs). The goal was to use the fractal properties of reservoir parameters to create a reservoir model with representative values of rock properties at scales below the resolution limit of the seismic technique (i.e. about 30 m) in order to discover whether including variability at the smaller scales would allow such models to perform better. One of the advantages of AFRMs is that their structures, heterogeneities and anisotropies can be accurately controlled. Al-Zainaldin et al. (2017) generalised to three dimensions a workflow and associated coding which had been developed for modelling fractures with rough surfaces in two dimensions, some time previously, and that had been available online for researchers as the software Synfrac (Brown et al. 2011;Glover et al. 1997Glover et al. , 1998aOgilvie et al. 2002Ogilvie et al. , 2003Ogilvie et al. , 2006. About the same time, other workers were also looking at 2D fractal objects for modelling transport properties within a reservoir (Le Ravalec et al. 2000). Glover et al. (2018) used the Al-Zainaldin et al. (2017) approach to compare the efficacy of AFRMs against conventional seismic data-mediated krigged reservoir models, and discovered that while there was no advantage for homogeneous, isotropic reservoirs, the AFRMs performed startlingly better for heterogeneous and anisotropic reservoirs.
The aim of this paper is to use generic AFRMs (GAFRMs) with well-defined heterogeneities and anisotropies to study the effects of well placement and, critically, injector-producer orientation on reservoir production parameters in reservoirs with significant heterogeneity and lateral anisotropy. Production from the GAFRMs has been simulated using a 6-spot vertical injector-producer well pattern as a function of orientation of the well pattern with respect to the reservoir, with simulations made every 20° of relative rotation. We have calculated average oil and water production rates (AOPR and AWPR, respectively), cumulative oil production (COP) and water cut as well as the time to water breakthrough using reservoir simulation. These parameters were calculated for every 20° of rotation of the resulting reservoir models with respect to a pre-defined and constant well pattern containing injector and producer wells.
It is important to understand that the modelling carried out in this paper is not specific to any given field. It is generic and has the advantage that the heterogeneity and anisotropies can be completely controlled. This means that we can examine in detail the effect of changing a single parameter. It also means that we do not suffer from missing data or need to tune the model to historical production. The parameters are either those we wish to vary in order to understand the model more clearly, or standard typical values that represent a typical reservoir. In the first case, an example would be the variation of horizontal anisotropy, while in the latter case an example would be the values of oil and water viscosity we have used. Figure 1 shows the generalised procedure followed for each set of simulation runs carried out in this work. The left-hand branch of the figure shows the creation of a fully specified AFRM with a pre-defined size, structure, heterogeneity and anisotropy defined by the input properties. The AFRM can be used directly in reservoir simulation. The right-hand branch of the figure considers the creation of a set of 18 well patterns created by the rotation of a six-spot injector-producer unit cell by increments of 20°. Simulation was then carried out to give the production rate and cumulative production of all reservoir fluids, as well as the reservoir fluid pressure as a function of spatial position in the model and of time.

Implementation of Advanced Fractal Reservoir Models
The study was conducted using Advanced Fractal Reservoir Models (AFRMs) (Al-Zainaldin et al. 2017;Glover et al. 2018). These models take full account of partial spatial correlation of the petrophysical parameters which compose the model as well as the fractal nature of the Corey exponent (λ), which defines water distribution and relative permeabilities. It is a property of AFRMs that they contain information at all scales, and hence are extremely well suited to describe heterogeneous and anisotropic reservoirs where variability at small scales often controls fluid flow at larger scales and vice versa. This work uses generic AFRMs in which the petrophysical parameters are normalised to realistic typical values commonly found in clastic reservoirs, but without conditioning the models to represent any given reservoir structurally, although the latter is now possible (Glover et al. 2018).
The AFRMs were created using modified Matlab codes that were developed by Al-Zainaldin et al. (2017). The modelling approach derives originally from the theory of Mandelbrot and Van Ness (1968) which described fractional Brownian Motion (fBM). In our case, we use a Fourier filtering method to simulate fBM. Initial codes for 2D rough surfaces (Saupe, 1988) were first modified in order to model 2D unmatched rough fractures (Brown et al. 2011) and then partially matched 2D rough fractures (Ogilvie et al. 2006;Glover et al. 1997;1998a, b). Al-Zainaldin et al. (2017) extended the fractal 2D method to three dimensions so it could be used as a base for generating 3D reservoir models. The resulting AFRM code is capable of creating three-dimensional fractal geological models of reservoirs with pre-defined means and spreads of common petrophysical parameters, including (i) porosity (ϕ), (ii), cementation exponent (m), (iii) characteristic grain size (d g ), (iv) the effective permeability (k eh = k ex = k ey and k ev = k ez ), (v) capillary pressure (P cap (S w )), (vi) water saturation (S w ), and (vii) relative permeabilities (k rh = k rx = k ry and k rv = k rz ), all with pre-defined heterogeneities, defined by a set of fractal dimensions, and horizontal and vertical anisotropies, defined by a set of anisotropy factors. Although the reservoir models are fully deterministic (i.e. can be repeated exactly from the same data set), they also include a random number seed which acts as a key which defines the individual structure of the reservoir, as fully described in Al-Zainaldin et al. (2017) and Glover et al. (2018).
The heterogeneity of a reservoir is given by the fractal dimensions which define the input petrophysical parameters. In our case these petrophysical parameters were (i) Fig. 1 The general modelling procedure used in this work. The left branch creates AFRMs with different structures of heterogeneities and anisotropies. The right-hand branch defines the position of all of the wells for each set of measurements. Simulation brings the two branches together to produce the production rate, cumulative production fluid pressure and water cut data as a function of spatial position in the model and time. Example input parameters: Seed = 73977373, = 3.1, χ yz = 4, χ xy = 0.8 porosity, (ii) cementation exponent, (iii) grain size, and (iv) the Corey exponent (λ), with permeability and all other reservoir model parameters being calculated from these input data. However, such modelling could be carried out with other input parameters, for example using a fractal volume to model permeability directly. It is important that the reservoir has some heterogeneity because a completely homogeneous reservoir cannot be anisotropic. Equally, a reservoir that is extremely heterogeneous approaches spatial distribution which is random, and such a random distribution in three dimensions is also isotropic. The moderate heterogeneities that are commonly found in hydrocarbon reservoirs can be significantly anisotropic.
Anisotropy is defined in the models by a set of three anisotropy factors χ xy , χ yz and χ zx , which have rotational symmetry in such a way that the definition of any two factors automatically defines the third. The anisotropy factor χ xy defines how anisotropic or stripey the reservoir is in the xy plane, with χ xy > 1 leading to stripes parallel to the x-direction, and χ xy < 1 leading to stripes parallel to the y-direction, such that if χ xy is allotted some value n providing a given degree of striping is in the x-direction, allotting χ xy the value 1/n provides the same degree of striping in the y-direction. The other two anisotropy factors define vertical anisotropy, which is generally associated with horizontal or subhorizontal layering in the reservoir. Commonly, we set χ yz > 1 to provide different degrees of vertical anisotropy giving striping parallel to the y-axis, but we could get the same result by setting χ zx < 1. Figure 2 shows a complete AFRM. The details of the construction of models such as these are given in Al-Zainaldin et al. (2017) and Glover et al. (2018). In this particular model, fractal volumes for porosity, grain size and cementation exponent were created with individual fractal dimensions to define the heterogeneity of the model ( = 3.1), and an imposed horizontal and vertical anisotropy of χ xy = 4 and χ yz = 0.8. The histograms for each of these conditioned fractal distributed volumes are also given, and can be seen to be approximately normal distributed with a mean value and standard deviation equal to that imposed by the modelling software for each parameter ( ̄ = 0.18. σ ϕ = 0.03, m = 1.7, σ m = 0.05, d = 1.2 × 10 −4 m, σ d = 0.55 × 10 −4 m). The calculated permeability map is log normally distributed with a modal value of about 75 mD, and the calculated water saturation map varies between about 0.16 and 0.76. Table 1 shows the parameters used in the AFRM modelling for all of the models simulated in this paper; where a single value is given, this has been held constant for all simulations. The mean and standard deviations of porosity, cementation factor and characteristic grain size were taken as typical values for a clastic reservoir, but would be obtained from core data in implementations modelling specific reservoirs.
It is important to note that in this work the same degree of horizontal (χ xy ) and vertical (χ yz ) anisotropy is implemented for the normalised fractal volumes that represent (i) porosity, (ii) cementation exponent, and (iii) characteristic grain size. The model allows individual control of the anisotropies of each, if required. The consequence of anisotropy in these initial three parameters is that there is anisotropy in the calculated fractal permeability volume through the application of the RGPZ equation Rashid et al. 2015b), in the calculated synthetic poroperm plots, calculations of capillary pressure, relative permeability curves and water saturations (Al-Zainaldin et al. 2017).
Recognising that there is also a microscopic component to permeability anisotropy, as discussed in full in the next subsection, the horizontal permeability of a given cell in all lateral directions was taken as equal to that given by the calculated fractal permeability volume, while the vertical permeability of a given cell was taken as equal to that given by the calculated fractal permeability volume divided by 5. This is the same as recognising a sub-cell scale k v /k h ratio equal to 0.2.

Anisotropy
This subsection considers anisotropy in more detail since this paper is concerned primarily with anisotropy and directional flow. Anisotropy occurs at different scales in real reservoirs and in the modelling that we implement. At small scales, the effect of mineral morphology, such as micas, or compaction, as in shales, causes there to be implicit anisotropy in the rock. Generally, such anisotropy leads to differences in fluid flow vertically and horizontally, with vertical permeability commonly being 5 to 10 times less than horizontal permeability. More often than not, this implicit anisotropy is not developed in the horizontal plane, with permeabilities not changing azimuthally. In some situations however, it is possible to observe small-scale azimuthal anisotropy in rocks which have been subjected to significant differential horizontal stresses. In this work, the small-scale implicit anisotropy in permeability has been set to k ix = k iy = 10k iz (where the subscripts refer to implicit anisotropy and the Cartesian direction, respectively) for all modelling, effectively removing small-scale implicit anisotropy from the analysis of our results.
At larger scales, explicit anisotropy is developed due to the juxtaposition of rocks with different textures even if they may be lithological and mineralogically identical. This larger scale anisotropy is defined by the explicit spatial arrangement of zones of different petrofacies in all three directions. Commonly, it is mainly controlled by bedding, at scales varying from below that at which the gamma ray log might measure to greater than the seismic Fig. 2 a The physical property map of fractal data volumes conditioned to represent porosity and permeability using heterogeneity ( = 3.1) and isotropic (χ xy = 4, χ yz = 0.8). b The physical property map of fractal data volumes conditioned to represent water saturation, cementation exponent and grain size, using heterogeneity ( = 3.1) and isotropic (χ xy = 4, χ yz = 0.8) resolution. However, lateral variations also have a major part to play in many depositional environments. A good example of this can be found in Al-Qassab et al. (2000), where lateral variation is as a result of variations in sabkha activity.
In simple layer cake modelling, explicit anisotropy in the z-direction arises from the harmonic mean of the value of k iz for each layer and taking into account the thickness of each layer. In conventional reservoir models, the discrete element numerical simulation will take account of the explicit anisotropy in every direction. However, the explicit anisotropy will not be quantified nor reported.
In this work, the explicit anisotropy is well defined mathematically in each of the Cartesian directions. The AFRMs used in this work required the three-dimensional inverse Fast Fourier Transform of a four-dimensional spectral density function (Al-Zainaldin, 2017). The slope of the spectral density function defines the heterogeneity, or fractal dimension, in each direction. The power function can also be manipulated mathematically so that the resulting AFRM exhibits anisotropic behaviour. Brown et al. (2011) presented a generalised spectral density function to generate anisotropic 2D fractal surfaces analogous to their isotropic form. Extending this to the three-dimensional case, we multiply each frequency component that corresponds to a certain spatial direction by a constant factor, which we define as the anisotropy factor χ that was mentioned earlier in this paper. This results in a generalised spectral density function where the amplitude spectrum coefficient α klw follows a power law behaviour defined by the parameter H, which is related to the fractal dimension by = 4-H, and where k, l and w is the position in the spectral density function (analogous to x, y, and z), and the three factors χ xy , etc. are the anisotropy factors mentioned earlier in the text.
The permeability for each cell in the three-dimensional matrix is calculated from the three-dimensional volume representing the original three normalised parameters (porosity, cementation exponent and grain size) using the RGPZ equation Rashid et al. 2015b). This parameter can be normalised by correlation with core measurements or normalised by the implicit anisotropy to give separate vertical and horizontal permeabilities. Consequently, we obtain a full and well-defined effective permeability tensor for each block in the model (Duquerroix et al. 1993;Noetinger and Haas 1996). Schneider et al. (2018) have shown very nicely how different discretisation schemes can lead to unexpected solutions that do not agree with each other. The advantage of the AFRM approaches is that the model is cell-centred, providing all implicit and explicit permeabilities for each cell. Variability from cell to cell is not sudden, as a random spatial distribution would be, but follows the natural variability that is associated with the imposed fractal heterogeneity implicit in the modelling.
Both the fractal dimension (heterogeneity) and the anisotropy values in the final normalised fractal volumes were confirmed by the use of independent analytical software. In this work the small-scale, implicit, anisotropy is fixed with k iv /k ih = k iz /k ix = 0.1 and k iH /k ih = k iy /k ix = 1, which arises from the statement made in the first paragraph of this subsection. The large-scale value of k v /k h arises from the fractal structure of each AFRM. We believe that the anisotropy factors previously defined represent a way of quantifying these large-scale anisotropies, but the flow results from simulation would be their best macroscopic measure. A more mathematical treatment of anisotropies may be found in Al-Zainaldin (2017)

AFRM Simulation
A large number of AFRM simulations have been carried out. In each case a six-spot (injector-producer) well pattern was used, and was rotated with respect to AFRMs. All simulations were carried out using the Roxar Tempest (ver. 7.0.4) black oil simulator. All of the models tested consisted of 6,291,456 voxels (i.e. 256 × 256 × 96 voxels) with a 12.8 km × 12.8 km reservoir extent and 96 m thickness, with a reservoir top at 2000 ft and an oil-water contact at 2096 ft. The initial fluid pressure was 2300 psi and the temperature was 120 °C. Simulations were carried out as a function of time for a maximum of 30 years, and with a maximum fluid production rate of 5 Mstb/well/day. In each case oil average or production rate, any gas production rate as well as the cumulative volumes of each fluid were calculated, together with the time to water breakthrough at the producer wells. These calculations were carried out for the rotation of the six-spot well pattern with respect to the AFRMs for rotational increments of 20°, noting that it is much easier to implement the rotation of the six-spot well pattern than it is to rotate the entire AFRM. Care was taken to ensure that the AFRM was significantly larger than the well spacing in order that variable edge effects in the results are not significant.
Table 1 also shows the seed values used in the modelling. While these are not necessary for a full understanding of the paper, they would be necessary for anyone wanting to repeat the modelling. Each seed contains enough information to create the original structure of the model being tested. In this work, some models were tested on just one structure for which the random number seed is given. However, it was recognised that there is no guarantee that this one structure would necessarily be representative of any given anisotropic reservoir. For example, by chance, it might contain repeating structure that would not be considered typical, or a predisposition towards either high or low permeabilities towards either the injector or the producer. Any of these chance situations might cause the results of the modelling to be, if not invalid, then irrelevant. Consequently, the modelling was also carried out on a set of seven reservoirs with different structures. The seed values for each of these are given such that the reservoir modelling can be repeated if required. In this modelling, only the seed value was changed, leaving the fractal dimension, anisotropies and all other parameters constant. Consequently, mean simulated values from the seven different structures were calculated and presented, together with an error bar indicating their spread. This modelling, based on a cohort of AFRMs, is much more likely to provide results which are typical of the behaviour of the reservoir.
It is recognised that, although simple to implement, the use of rectangular grids may induce apparent anisotropy in the modelling results as an artefact (Ababou 1996). Such anisotropy can sometimes be seen in results from uniform rectangular grids populated with homogeneous and isotropic reservoir properties, but is usually small. In this work, we first tested the AFRM approach with homogeneous and isotropic AFRMs and observed no significant grid-induced anisotropy. The very small variations that were observed were slight preferences to fluid flow in the injector-producer direction. However, these were 100 times smaller than preferential flow induced by even the smallest heterogeneities and anisotropies used in this work. Referring to Table 1, the water-to-oil viscosity ratio was 0.001752, while the mobility ratio was calculated to be over 200, which would make the displacement naturally unstable and amplify any anisotropy affect. Consequently, the results of this paper should be seen as deliberately the worst-case scenario.

Modelling Results
The results are presented in three parts. In the first part the directional behaviour of a heterogeneous but isotropic reservoir will be presented and discussed. This reservoir has a fractal dimension of = 3.1 and anisotropy factors χ xy = χ yz = 1. In the second part, the directional behaviour of a heterogeneous reservoir with different degrees of anisotropy will be given and analysed. Both of these parts involve simulation with a single reservoir structure defined by a single random number key value. This reservoir also has a fractal dimension of = 3.1 and variable anisotropy factors χ xy = 1, 0.9, 0.8, 0.7, 0.6 and 0.5, and χ yz = 4 or 5. Since there is the potential for variability in the results to depend on different reservoir structures, we recognise the necessity to carry out simulations for a number of different structures with different random number key values. Consequently, the final part will consider the mean behaviour of simulations made on reservoirs as a function of direction and anisotropy, as well as noting the maximum and minimum simulated values such as the envelope in which reliable simulated parameters are likely to fall.

Isotropic Heterogeneous Reservoir
This simulation used a single reservoir structure with a fractal dimension of = 3.1 and anisotropy factors χ xy = χ yz = 1. The purpose of the simulation was to analyse the production behaviour from an isotropic heterogeneous reservoir in order to confirm the expected lack of directional behaviour in the simulated outputs. Figure 3 shows the simulated oil production rate as a function of time starting at an arbitrary year, namely 2000. The six curves each show the production for the scenario where the six-point injector-producer pattern is rotated by a given amount with respect to the reservoir. It might be thought that there would be an insignificant difference in the time profiles for the production rate because the structure is deliberately isotropic. However, this is not the case. The separation in the curves, particularly in the first half of production indicates that there is, for this reservoir, directional behaviour. However, close analysis of the order of the curves with respect to their rotation shows that it is not a case of there being a simple preferred direction, but that there is significant variability associated with any rotation.
We hypothesise that the apparent directional behaviour is due to the heterogeneity of the reservoir. The act of rotating the injector-producer pattern with respect to the reservoir causes injectors and producers to be placed in new positions which may be better reservoir quality or worse reservoir quality than the previous rotational position. Consequently, and without any imposed anisotropy to the reservoir, one would expect variability in the simulated production values. Since, for each rotation, the overall effect of whether each of the injectors and producers is in better or worse quality reservoir is not known or controlled, it is expected that there would be a variability in the production parameters for each rotational position which is defined by heterogeneity rather than by anisotropy.
It should be noted that the behaviour in Fig. 4 occurs in three parts. In the first part (approximately years [2000][2001][2002][2003][2004], the reservoir is producing primarily due to its initial formation pressure and natural drive which decreases from initially high values, and in our simulations the oil production rate was capped at 15 Mstb/day. From 2004 to about 2030, the decrease in production rate is allayed and even increased by the injection of water at the injector wells which maintains reservoir fluid pressure. In later production, all curves die away approximately exponentially as the reservoir becomes depleted. It is the injectionsupported production which is most affected by the rotation of the injector-producer pattern with respect to the reservoir structure, which is consistent with the hypothesis set out earlier.

Anisotropic Heterogeneous Reservoirs
The results presented in this section are all the results of simulations of synthetic reservoirs which have a fractal dimension of = 3.1, as before, but with the addition of variable horizontal anisotropy factors χ xy = 1, 0.9, 0.8, 0.7 and 0.6, and vertical anisotropy factors χ yz = 4 or 5.
The imposition of an externally defined anisotropy within the reservoir would be expected to lead to anisotropic flow, with either preferential oil or water flow at different stages in the reservoir production. Simulations have confirmed this. We first consider the effect of varying the vertical anisotropy and then consider the effect of varying the horizontal anisotropy. Fig. 3 a Six-spot well arrangement with three injector and three producer wells. b The rotation of the wells with respect to any AFRM under test. c Schematic diagram of the method for rotating any point (x,y) around any given origin to a new position (x',y') through an angle of ( − ) The tests carried out in this paper all considered six vertical injection production wells rotated about a vertical axis. It would be expected that anisotropy in the vertical direction leading to apparent layering of reservoir quality defining parameters would have little effect on production parameters. We have carried out simulation tests for a range of vertical anisotropy and found that varying this parameter has negligible effect on the hydrocarbon production rate and cumulative production as well as the water cut and breakthrough times. There are several reasons for this. First, we use a strongly anisotropic discretisation mesh between the horizontal and vertical directions, which gives rise to large differences between vertical and horizontal transmissibilities. However, it is also important to consider that all of our tests were carried out with vertical injector and producer wells that were perforated throughout the modelled volume. Consequently, there is little scope for vertical flow unless channelled by local anisotropy. The large vertical anisotropy factors we use in all of our reported modelling, χ yz = 4, ensures that the reservoir is significantly layered, and that vertical flow is not important such that we can concentrate on the effect of azimuthal anisotropy. We would expect much more of an effect if a set of horizontal wells were being used. In such a case, the chance positioning of one or more horizontal wells in either good or bad quality reservoir would make a significant change in production parameters obtained by simulation.
Accepting that the horizontal anisotropy will exert a main control on horizontal fluid flow between vertical wells, most simulation was carried out for a heterogeneous reservoir with a fixed given vertical anisotropy and a variable number of different horizontal anisotropies (i.e. χ xy = 1, 0.9, 0.8, 0.7 0.6 and 0.5). Fig. 4 Oil production rate as a function of time for the full simulation of a single heterogeneous ( = 3.1) isotropic (χ xy = χ yz = 1) reservoir when produced with a six-spot injector-producer pattern rotated with respect to the reservoir structure between 0 and 90° in six steps Figure 5 shows the directional distribution of production rate 3, 10, and either 25 or 30 years after the beginning of production for a heterogeneous reservoir with a fractal dimension of = 3.1 and two different horizontal anisotropies. In both cases, the same random number seed keys were used (see Table 1), implying that the underlying structure of all the AFRMs is the same, with only the degree of anisotropy in the fractal volumes changing. The two parts of the figure on the left-hand side (Fig. 5a, d) are for a layered horizontally isotropic reservoir (χ xy = 1, χ yz = 4), while those in the middle (Fig. 5b, e) are for a layered reservoir with a small horizontal anisotropy (χ xy = 0.8, χ yz = 4) and those on the right (Fig. 5c, f) are for a layered reservoir with a larger horizontal anisotropy (χ xy = 0.5, χ yz = 4). The rose diagrams show changes in the directional production as a function of time and clearly illustrate that a reservoir controlled by less heterogeneity results in higher and longer production (i.e. to 2030).
Taking the horizontally isotropic case first (Fig. 5a, e), it is clear that the heterogeneity of the reservoir is inducing some directional behaviour aligned with the 100-280° axis, as well as a preferential flow direction towards about 160°. Initial production is preferential in these directions. These directions occur by chance and arise due to the interaction between the injector-producer well pattern and heterogeneity of the rock. Consequently, as in the first results in this paper (for the homogeneous |AFRMs), directional behaviour can be induced by heterogeneity without it necessarily being due to anisotropy in the reservoir parameters. Such directional behaviour is not predictable. It is worthwhile noting that the anisotropy in production rate decreases over time, with the late production rate data showing not only an overall smaller production rate, as expected, but also demonstrating less directional sensitivity. However, vestiges of the original directional sensitivities still remain even in the case of the late data. Taking the slightly anisotropic case (Fig. 5b, e), we note that the horizontal anisotropy factor used here (χ xy = 0.8) induces only a mild anisotropy in the reservoir in a direction parallel to the y-axis. The preferential flow is clearly aligned with the 90-270° axis as expected from the imposed anisotropy, but since the heterogeneity of the reservoir itself prefers this direction it is difficult to tell the extent to which the imposed anisotropy is forcing fluid flow. This is a difficulty associated with doing measurements on a single reservoir, and is addressed later in this paper by using a cohort of reservoirs. However, we present the data because we feel that, taken with other simulations, it is instructive. The preferential flow lobe at 160° in the isotropic case is less pronounced in the production rate for the slightly anisotropic case, and we hypothesise that the anisotropy is beginning to attenuate this directional flow artefact which is caused by heterogeneity alone.
The case with a high degree of anisotropy (χ xy = 0.5, Fig. 5c, f) represents an extension to the previously perceived behaviour. Once again, production rate for all years is smaller than for the other two cases. Once again initial production is highly anisotropic arranged in the 90-270° direction. The degree of anisotropy in the flow is greater than in previous cases and the flow in the 160° direction has been completely overwhelmed by the imposed anisotropy of the flow system, with a slight bump in the curve the 2003 at 160° being all that remains. This highly anisotropic flow is still present by 2010, though at a lower production rate, as expected. By 2025 the flow shows vestiges of direction behaviour, but is much more isotropic as the reservoir approaches its end.
Consequently, Fig. 5 shows that imposed reservoir model anisotropy can have a very significant effect on flow anisotropy and that anisotropy can be effective over a significant fraction of the lifetime of the reservoir, and is consistent with the results of Al-Zainaldin et al. (2017). However, built-in directional sensitivities which are purely due to heterogeneity of the reservoir can also be significant in reservoirs which are isotropic or have a small degree of anisotropy, but will be negligible in reservoirs which have a large degree of anisotropy. Figure 6 shows the production rate as a function of direction and horizontal anisotropy. It is clear that whatever directional preference the isotropic reservoir exhibits, increasing the anisotropy in the 90-270° direction amplifies the effect, and in doing so increases the production rate. The reason for the flow in the 160° direction being increased as anisotropy increases is currently unknown, but may arise as an artefact of carrying out simulations on a single reservoir in this case. Recognising this problem, simulations have been carried out on a cohort of reservoirs sharing the same physical properties but having different random number seeds and therefore different structures in the last part of this paper.

Water Production
The water cut as a function of direction and horizontal anisotropy at two different dates (20 years and 30 years after production starts) is given in Fig. 7. Looking at the 2020 data, significant anisotropy in a water cut is apparent for all values of horizontal anisotropy in the reservoir properties, including that belonging to the isotropic reservoir (χ xy = 1). As horizontal anisotropy increases (i.e. lower values of χ xy ), anisotropy in the water cut also increases, leading to highest values of water cut occurring in the direction of anisotropy (90-270°), reaching 25% in 2020 for the highest degree of anisotropy, which occurred at 280°, compared to a value of 18% for the isotropic reservoir, which occurred at 80°. This result concords well with experience in heterogeneous reservoirs, where water flow anisotropy occurs as a result of preferential flow channels forming between injector and producer wells, leading to larger water cuts at the producing wells and much earlier times to water breakthrough.
By 2030, all water cut values are much higher, as would be expected, but with only traces of the anisotropy remaining. For late-stage production, the water cut is significantly higher for anisotropic reservoirs than the isotropic case, suggesting that replacement of well-connected hydrocarbons with water in early and middle production has developed water flow pathways which do not have a preferred direction in the reservoir despite the reservoir having anisotropic distributions of porosity, grain size, cementation exponent and permeability. Fig. 6 Simulated production rate of an AFRM with the same random number seed structure key, heterogeneity ( = 3.1), vertical anisotropy factor (χ yz = 4) but varying the horizontal anisotropy factor (χ yz = 1, 0.8, 0.6) In this work, we observe that the preferential flow of injected water through highly permeable regions causes the deviation of the waterfront towards more highly permeable zones, increasing the production of oil from those zones significantly but by-passing oil in less permeable zones. The results from our 3D AFRM modelling approach is compatible with the experimental study presented by both Ahmedi et al. (2019) and Ali and Al-Qassab (2000) in this respect.
A small amount of modelling work with AFRMs models has also been done to examine this effect (Glover et al. 2018), where it is clear that the optimal production of oil occurs when both injector and producer wells are situated in high permeability rock.

Mean Behaviour of Heterogeneous Anisotropic Reservoirs
All the results discussed so far in the paper have been obtained using a single random number key value, which results in a single specific reservoir structure. There is, of course, no guarantee that such a structure does not contain a particularly unusual structure which by chance has a particular preferred flow direction other than that imposed by the modelled anisotropy factors. The flow in Fig. 5 towards 160° was recognised as possibly arising from such behaviour. Fig. 7 Results from reservoir simulation of water cut for isotropic and anisotropic reservoirs with the same random number seed structure key, heterogeneity ( = 3.1) and vertical anisotropy (χ yz = 4) but 4 different horizontal anisotropies (χ xy = 1, 0.8, 0.6). a Data for 2020 (20 years after production start), and b for 2030 We have carried out modelling and simulation for a number of different reservoirs sharing the same physical properties (dimensions, heterogeneity, anisotropy factors, mean and spread of modelled parameters of porosity, cementation exponent and grain size, etc.) and differing only in the random number key that defines the specific reservoir structure. Once again, simulation has been carried out as a function of time and 18 rotations of the injector-producer well pattern with respect to each reservoir. Such modelling is computationally extremely intensive. Consequently, a cohort of only 7 reservoirs was possible to be considered. Nevertheless, the mean and bounding behaviours of the production parameters from the results of simulating these 7 reservoirs are sufficient to define a reliable envelope of simulated production parameters and to exclude chance reservoir structures that might cause out-liers in the simulated results. Figure 8 shows the mean simulated production rate data from simulating the 7 reservoirs as a function of direction and time together with the maximum and minimum behaviours which are shown by the error bars. The results clearly show that using reservoirs with different specific structures (defined by different random number keys) has an effect on the simulated parameters, as would be expected. However, the variability resulting from using different structures is surprisingly small as indicated by the size of the error bars, and is small compared to the variability encountered as a function of rotation or of time. This result gives us confidence that we can use simulation on a single reservoir to give a reasonable estimation of the production behaviour. Figure 8 presents an ellipsoidal shape with the long axis arranged in the 80-260° direction, which accords well with the imposed horizontal anisotropy in the 90-270° direction. Application of Rayleigh's test shows that the data are anisotropic and that there is less than a 2.21% probability that such an alignment would occur by chance. Furthermore, Fig. 8 shows no particular lobes indicating preferential flow in other directions as was the case in Fig. 5. Both the early and mid-term productions exhibit anisotropic production rate in the expected direction, with the early production rate being larger than mid-term production rate, again as expected. The late production rate is lower yet again, but has become isotropically distributed in accordance with the earlier results and the expectation that reservoir Fig. 8 Mean production rate together with maximum and minimum values in production rate for 7 reservoirs with different random number structural keys but the same physical parameters, dimensions, heterogeneity ( = 3.1), vertical anisotropy factor (χ yz = 4) and horizontal anisotropy factor (χ xy = 0.6) as a function of direction at different dates during production anisotropy does not control reservoir behaviour when all of the initially produced hydrocarbon has been replaced with water, leaving once highly connected hydrocarbon paths now occupied by water and hence hindering late-stage hydrocarbon flow.

Discussion
The majority of the results shown in the section above are intuitive and agree with what we understand from experience. For example, (i) that predominantly horizontal flow between the vertical wells is not greatly affected by varying vertical anisotropy, (ii) that production rate and cumulative production are both strongly affected by orientation to horizontal anisotropy as well as heterogeneity, and (iii) we would expect that heterogeneity would influence the results depending on the positions of the well to some extent even if no anisotropy is present.
The strength of this paper is not that it presents material which disagrees with experience or field observations, but that this new modelling approach provides valid and mathematically rigorous modelling, which is entirely consistent with past observations, and is also generic (i.e. it does not rely on making a value judgement on a small number of individual reservoirs or cases that we label as experience). Moreover, the parameters describing heterogeneity and anisotropy which have been developed and used in this AFRM modelling are in themselves mathematically rigorous, which has not previously been the case, and may be varied in future generic AFRM modelling to examine a range of heterogeneities and anisotropies present in a particular reservoir. Hence, the paper can also be viewed as a validation of the AFRM reservoir modelling approach using existing understanding of anisotropic reservoirs as a gold standard.
However, such generic modelling also has the potential for discovering associations and processes which have not been previous realised, or are otherwise counter-intuitive. An example of this is that the optimal production from isotropic heterogeneous reservoirs was found to occur when both injection and production wells are sunk in high reservoir quality rock, where one would expect intuitively that at least the injection wells, and perhaps the production well, should be sunk in moderate reservoir quality rock in order to improve the sweep efficiency of the injection process and reduce the volume of by-passed oil. Comparing the production rate data in Figs. 5 and 6, which are for a specific position of injector and producer wells in a single AFRM with that of Fig. 8, which is for the mean of 7 structures, indicates that the effect of heterogeneity can have an effect at least as large as welldeveloped anisotropy.
This effect has been confirmed by the results of an associated study which examined the placement of injector and producer wells in near-homogeneous ( = 3.1) and heterogeneous ( = 3.9) but otherwise isotropic AFRMs. In this set of reservoir modelling it was the placement of injector and producer wells that was studied. Tests for well placement are difficult to make because of the large number of possible combinations of scenarios. A total of 810 simulations were carried out testing the effects of combinations of the following scenarios, where the total number of possible combinations is 1350: • Distribution of both the injector and producer wells randomly using 3 random configurations. • Distribution of the injector wells randomly, but deliberately placing the producer wells in high permeability rock with the exact location allocated at random.
• Distribution of the producer wells randomly, but deliberately placing the injector wells in high permeability rock, with the exact location allocated at random. • Deliberate placement of both the producer and injector wells in high permeability rock, with their exact position. • Repeating the last three scenarios above, but with any deliberate placement being in low permeability rocks. • All of the tests above for combinations of 1 to 10 producers and 1 to 5 injectors. • Variable heterogeneity represented by = 3.1, 3.5 and 3.9.
In all cases the placement of wells was constrained such that wells had to be further apart than a minimum inter-well distance (d min ). Discrimination between high and low permeability rock was taken to be k cut-off = 120 mD. Fig. 9 Average oil production rate (AOPR) for combinations of 1 to 10 production wells and 5 wells for near-homogeneous ( = 3.1) and heterogeneous ( = 3.9) AFRMs. In each case, tests are carried out for 3 random well patterns of injectors and producers, with producers placed deliberately in low or high permeability regions of the model, and with both injectors and producers placed in high permeability zones. The distinction between low and high permeability for the purposes of this test was set at k cut-off = 120 mD Figure 9 shows some of the results for the scenario with 5 injector wells for the nearhomogeneous and heterogeneous case ( = 3.1 and 3.9). The variation in data shown by the three random placements set the general variability that might be expected from random well placement. For the near-homogeneous reservoir, the difference between the 3 random cases is large for one production well, and becomes smaller as the number of production wells increases. This is consistent with the idea that additional randomly placed wells will sample the AFRM more representatively. For the near-homogeneous reservoir, the placement of production wells leads to insignificant improvement or degradation in a production rate compared to the 3 random reference cases as the number of production wells increases from 1 to 10. The best production rate is obtained when both injectors and producers are placed in high permeability rock, irrespective of the number of production wells used, but the advantage is not statistically significant.
For the heterogeneous reservoir, all AOPRs are significantly less than for the homogeneous case, as would be expected. However, the span between the 3 random cases follows the same general behaviour, large for one production well, and becoming smaller as the number of production wells increases, for the same reason as before. The improvement in AOPR arising from the deliberate placement of production wells in high permeability rocks is significant from 1 to 4 producer wells, but the advantage is slowly lost, while deliberate placement of production wells in low permeability rocks seems to have no insignificant effect on AOPR, irrespective of the number of wells. It is clear, however, that when both injection wells and production wells are placed in high permeability rocks, there is a large significant increase in AOPR that is maintained as the number of producer wells is increased.
Recalling that this last set of tests is on isotropic AFRMs, it is clear that well placement, whether by chance or deliberately, may have significant effects on the reservoir production parameters.

Conclusions
The sensitivity of production parameters such as hydrocarbon production rate, cumulative fluid production and water cut to the anisotropy of heterogeneous reservoirs have been studied by creating generic advanced fractal reservoir models (GAFRMs) and simulating hydrocarbon production from them as a 6-spot injector-producer pattern that has been rotated with respect to the reservoir structure. Modelling and simulation has been carried out for single reservoir structures and for mean and bounding behaviours of a cohort of different reservoir structures as a function of time and direction as well as vertical and horizontal anisotropy.
The outcome results of 3D AFRM modelling concurs broadly with experience from field observations, intuition and previous theoretical and experimental studies as shown in Ahmedi et al. (2019). In a sense, therefore, some of the results are not new. However, the paper has shown not only that AFRM modelling is consistent with the previous science and our expectations, the explicit and full control over the quantitative modelling parameters representing heterogeneity and anisotropy that it contains allow these parameters to be compared between different reservoirs instead of being simply labelled as 'high' or 'low' heterogeneity or anisotropy as has previously been the case.
The main findings of the study are as follows: • Production rate and cumulative production of both hydrocarbons and water as well as the water cut are insensitive to both the degree and direction of vertical anisotropy when producing from vertical wells as in this study. • Production rate and cumulative production of both hydrocarbons and water as well as the water cut are extremely sensitive to both the degree and direction of horizontal anisotropy when producing from vertical wells. • Defining the injector-producer direction parallel to horizontal anisotropy leads to high initial production rate, early water breakthrough and lower long-term cumulative hydrocarbon production compared to defining the injector-producer direction perpendicular to horizontal anisotropy. • Reservoir horizontal anisotropy causes anisotropy in hydrocarbon production rate and water cut in early and middle production, with the directional distribution of these parameters becoming isotropic in late production. • Early water breakthrough occurs when the injector-producer direction is parallel with the horizontal anisotropy of the reservoir, as injected water migrates much more easily along the high permeability channels. • Apparent directional behaviour can occur in heterogeneous reservoirs which do not exhibit explicit anisotropic reservoir properties because the location of the injectors plays a very important role in driving oil to the wellbore. Consequently, chance placement of injectors and producers in high permeability rocks can have the same effect as anisotropy even though no anisotropy is present.
The modelling carried out in this paper is of generic nature, and therefore not applicable to any given reservoir directly, but pertaining in general to reservoir behaviour and giving us information about how anisotropy affects flow. However, we have recently made progress in creating the mathematical basis and coding necessary for the conditioning of AFRMs to represent real reservoirs, which will be published in the future.