Efficient History Matching of Thermally Induced Fractures Using Coupled Geomechanics and Reservoir Simulation

Abstract

Waterflooding is a common recovery method used to maintain reservoir pressure and improve reservoir oil sweep efficiency. However, injecting cold water into a reservoir alters the state of in-situ formation stress and can result in the formation fracturing. In other words, it can cause the initiation and growth of thermally induced fractures (TIFs), even when the original fracture propagation pressure is not exceeded. TIFs can cause non-uniform distribution of the fluid flow in wellbores, a reduction in sweep efficiency, and early water breakthrough in nearby production wells. Modelling and history matching workflows that consider the dynamic nature of the TIF problem are critical. These workflows improve and validate reservoir and geomechanical models, identify and confirm observed TIF onset and propagation periods, and provide a history-matched sector model with the rock mechanical and thermal properties and stress gradients that can be used with confidence for subsequent studies. Modelling and the underlining assumptions of fluid flow in the TIF and reservoir matrix, as well as geomechanical changes due to cooling of the reservoir during injection, are detailed below. A 3D reservoir simulator coupled with 2D finite element TIF and geomechanical models were used to manually history match an injector (NI6) in the N Field sector reservoir model in which a TIF was observed. In this study, history matching workflows were developed to consider the dynamic nature of TIF development during waterflooding. The reservoir and geomechanical models were improved and validated via the observed TIF onset and propagation periods. The history-matched models produced can be used with confidence in subsequent studies. The practical workflows and guidelines developed here can be used in waterflooding operations during the modelling, design, and planning stages. The novelty of this study is the coupling approach of different complex processes done in order to capture dynamic changes during waterflooding operations. A similar history matching study could not be found in the literature.

1. Introduction

It is well established in the oil industry that a large proportion of water injection wells are fractured at some time in their extended life [1,2]. This includes injection wells initially injecting at a bottom hole pressure (BHP) lower than the original minimum horizontal stress (σ_hmin) [3]. When a cold fluid is injected into a hot reservoir, the rocks contract because the injected fluid’s temperature is significantly cooler than the reservoir temperature. This contrast in temperatures causes the in-situ stress to be reduced considerably. Numerical illustrations of the temperature and stress evolution around the injection wellbore are shown Figure 1 and Figure 2 for a square-shaped reservoir that includes a horizontal producer and vertical water injector [4]. Figure 2 illustrates that the minimum horizontal stress (σ_hmin) is reduced around the injection well due to the reduction in temperature. When the minimum horizontal stress (σ_hmin) falls below the BHP due to temperature changes, fractures may initiate and/or propagate [3,5,6]. This applies to both vertical and horizontal wells [2]. Fractures initiation and/or propagation resulted from thermal processes is referred as Thermally Induced Fractures (TIFs).

A real-world example illustrating these thermal phenomena is shown in Figure 3 [7]. The fracture gradient reduction is directly proportional to the higher temperature difference between the formation and injected water temperature in the Prudhoe Bay field for a number of step-rate tests. This shows the significant impact of a temperature difference between the formation and injected fluid on stress reduction, and hence on TIF initiation and/or propagation

Horizontal and vertical injection wells often face the possibility of creating TIFs. The following can result from TIF development in injection wells in waterflooding and non-waterflooding applications [8,9,10,11,12,13,14,15].

• Non-uniform outflow of the injected water from the well can move into the reservoir. The literature has showed that most of the injected water can be accepted by a small zone. This causes poor sweep efficiency and non-uniform pressure support.
• TIFs can grow out of a net pay zone that causes the injected water to be lost to the formation above or below the target zone. This results in an inefficient waterflooding processes and reduces oil recovery.
• TIFs can have a positive impact on disposal wells by improving injectivity and achieving injection requirements. However, TIFs propagating into caprocks and contamination of fresh water remain challenges.
• Sand production may take place earlier than predicted due to cooling reducing the temperature and thus rock strength. Sand production may be observed during back-production of injection wells during well shut-ins.
• TIFs can have a positive impact on water treatment facilities. Water quality specifications can be relaxed because of the influence of TIFs.
• TIFs can significantly improve injectivity in disposal and geologic CO₂ sequestration applications, allowing for the achievement of economic injection rates. Avoiding out-of-zone TIFs and the contamination of fresh water remain challenges.

Modelling TIFs can be very difficult, due to the complexity of the interactions among fluid flow in the fracture, geomechanical changes, and the reservoir matrix in real-world situations that require more advanced numerical models. The simultaneous coupling of fracture mechanics with geomechanical models and reservoir simulations is the only means of modelling these complex phenomena in an accurate manner. Modelling TIFs is a critical step in history matching TIF models for waterflooding projects that may experience TIFs. History matching is not common in the literature since the focus is normally on reservoir parameters (e.g., porosity, directional permeability, NTG, transmissibility, etc.) rather than the rock mechanical, thermal, and stress parameters. In this research, history matching workflows that consider the dynamic nature of the TIF problem have been developed. The ultimate objectives were to improve and validate the reservoir and geomechanical models, identify and confirm the TIF onset and propagation periods observed, and provide a history matched sector model with the rock mechanical and thermal properties and stress gradients that can be used with confidence in subsequent studies

1.1. History Matching for TIFs

The primary objectives of history matching include improving and validating reservoir simulation models and understanding the different processes occurring in reservoirs. History matching is generally performed in two stages [16]. The objective of the first stage is to match the average reservoir pressure and the second stage objective is to match individual well histories. Figure 4 summarizes the general strategy used for history matching. Although each reservoir is unique, these guidelines provide a first step in screening most reservoirs [16].

History matching involves adjusting a reservoir model until it reproduces the historical behavior of the reservoir. The variations between simulated and observed values in a reservoir are analyzed and the model parameters changed to obtain a good match. History matching is very common in the oil and gas industry and is used to obtain reasonable reservoir models for production forecasts.

There are two approaches normally used in history matching. The first is the manual method, which has the following steps [17]:

• Running simulations for the historical period;
• Comparing the results to actual field data;
• Adjusting the simulation input to improve the match; and
• Further analyzing and verifying the input data based on knowledge and experience.

The second approach is the automatic method, which involves the following additional features. It:

• Minimizes the objective function (i.e., the difference between the observed reservoir performance and simulation results and
• Excludes the human knowledge/experience factor, thus, the results can be unrealistic.

A different history matching procedure is investigated in this research. The dynamic nature of TIFs such that they change with shifting reservoir conditions creates a unique history matching problem. This problem requires the use of different processes, including geomechanics, fracture mechanics, thermal changes, and reservoir conditions, due to the need to consider the injection of cold water. Static reservoir parameters (e.g., porosity, directional permeability, NTG, transmissibility, etc.) are kept constant in this research, as they are in this case better understood and measured. It is the rock mechanical parameters, thermal properties, and stress gradients that are modified to obtain the final match, since a fracture mechanics problem is what is being investigated. Modifying the above parameters changes the TIF dimensions. This results in changes to the TIF’s transmissibility, which is altered based on the reservoir conditions until a good data match is obtained.

In this study, the manual approach was used to obtain the final match to the sector model of the N field provided by Operator E.

1.2. Previously Reported History Matching Studies with Thermally Induced Fracturing

Detienne et al [5] employed history matching for an injection well with a TIF in order to validate their analytical TIF model. The authors used a 2D hydraulic fracturing model to predict the TIF’s length and width. Thermo-elastic and poro-elastic effects were considered in the model, using the solution given by Perkins and Gonzalez [18].

Detienne et al. [5] described an offshore field in West Africa with 10 injection wells that had been injected with water for a period of three to five years. The wellhead pressure (WHP) remained fairly constant at 100–120 bars during this time. It was observed that the TIF resulted in an increase in the injectivity index by a factor of 1.5–2. In three wells, the injectivity index increased by a factor of 10, indicating a more substantial effect from TIF’s development. The WHP history matching for one of these three wells (Well A) via a radial flow model showed a good match during the first 30 days of Period A (see Figure 5). The injection rate was constant at approximately 200 m³/d.

The injection rate suddenly began to increase from 200 m³/d after 30 days, reaching 2000 m³/d after 120 days for Well A. Initiating a modelled TIF that was vertically confined to the height of the reservoir was found to give a good history match for the period between 30 and 60 days, as shown during Period B in Figure 6. However, the injectivity again increased suddenly after 60 days. It was found to be impossible to match the well history after 60 days if a reasonable set of reservoir parameters was employed, as shown during Period C in Figure 8.

The TIF height and reservoir permeability had to be increased in order to match the injection history after 60 days. The height of the TIF was increased from 20 to 120 m and later reduced back to 80 m (see Figure 7). The complete injection history could then be matched by the model after Point C (see Figure 8). The history matching workflow used in the case study faced the following disadvantages:

• The combination of varying TIF heights and reservoir permeabilities made it possible to match almost any potential change in flow rate and WHP, without proper physics checks.
• The 2D model used was not dynamic (i.e., the matching had to be performed at different steps for different regimes).
• The history matching was performed in an uncoupled manner.
• Changes to the geomechanical and rock mechanics were not considered.
• Future TIF propagation and growth could not be accurately predicted.

The above case study shows that a 3D reservoir model coupled with geomechanical and TIF models is necessary to simulate dynamic TIFs and overcome the above simplifications. Changes to the rock mechanical and thermal property values, as well as the in-situ stresses, all need to be considered to accurately model a TIF’s dimensions.

2. Problem Statement and Solution Approach

Almarri et al. [19] investigated TIF development, characterization, and identification in the real-world field cases utilized in the present work (i.e., N Field) using combined analytical and semi-analytical methods. The sector reservoir model received from the operator was already history matched by the operator with respect to the static reservoir parameters. The objectives of history matching in the present work were to:

• Develop workflows that considered the dynamic nature of the TIF problem,
• Improve and validate the reservoir and geomechanical models,
• Identify and confirm the observed TIF onset and propagation periods, and
• Provide a history-matched sector model with the rock mechanical and thermal properties and stress gradients that can be used with confidence in subsequent studies.

This study utilized a specialized reservoir simulator (i.e., REVEAL i.e., a commercial reservoir modeling software developed by Petroleum Experts) that enabled the integration of the reservoir, production, rock mechanics, and geomechanics processes in a single platform [4]. This simulator used a 3D reservoir simulation coupled with a 2D finite element (FE) TIF and geomechanical models to manually history match Injector NI6 in the N Field sector reservoir model where TIF was observed. This was performed by directly linking the numerical FE model for fracture initiation and propagation to the finite difference (FD) 3D simulation grids through connection factors. The coupling approach to reservoir, geomechanics, and fracture mechanics to history matching a thermally induced fractures is novel—a similar history matching study could not be found in the literature.

3. Methodology

History matching in reservoir simulation can be a time-consuming process that involves many trails. Including a TIF model adds complexity to this process. A workflow has been developed not only to obtain smooth and reasonable history matches for TIFs, but also to ensure that each step is correctly implemented. The workflow was as follows:

• The modelling equations and assumptions used in the history matching have been detailed to ensure that the essential elements of coupling among the fluid flow, heat flow, geomechanics, and fracture mechanics have been incorporated.
• The received model was converted from an isothermal to a thermal model. The pressure-volume-temperature (PVT) properties were defined at a single temperature. It was necessary to define the fluid PVT properties across the full range of temperatures encountered during the injection of the cold water into the reservoir.
• A simulation was run without TIF modelling and compared to the historical data. This step allowed for an examination of the performance of the injector and identification of any potential indication of a TIF. An increase in the injection rate corresponding to a reduction in BHP can be used as an indicator of a TIF formation, especially if there is a reduction in temperature in the fluid injected. The values calculated for the reservoir pressure around the injectors and bottom hole injection pressure did not match to reliable data at this stage.
• In-situ stresses (in this case, for the N Field) were estimated from the data obtained from offset wells supplied by the operator. This included an estimation of vertical stress, orientation of the minimum horizontal stress, and estimation of the magnitudes of the minimum and maximum horizontal stresses.
• Dynamic mechanical properties of the rock (such as the Young’s modulus, Poisson’s ratio, toughness, and Biot’s coefficient) were determined using either compressional and shear velocities from acoustic log data provided by the operator, or when necessary, by employing published correlations.
• A TIF was then modelled to assess uncertainties and data reliability, as well as to investigate the possibility of TIF initiation and propagation using parameters obtained in steps iv and v. The results were then compared to the historical data.
• The history matching was then finalized after rectifying uncertain parameters (e.g., reservoir pressure and BHP). Workflows at different levels of the history matching were then used.

4. Thermally Induced Fracture Modelling

4.1. Introduction and Underlying Assumptions

One common method of modeling hydraulic fractures within a reservoir simulator is to modify the flow in an FD reservoir by changing the grid properties to thin blocks with high porosity and permeability in order to simulate fracture flow. Alternative methods adopted for hydraulic fractures have included local grid refinement and equivalent effective wellbore radius methods. These are not suitable for TIF modelling for the following reasons [20]:

• The importance of flow in the reservoir for long-term TIF growth modelling is not appreciated.
• Large time steps are required for TIF modelling.
• TIFs have much higher leak-off rates.
• Thermal aspects are normally neglected in hydraulic fracture models (i.e., the thermo-elastic stress is not considered).

Therefore, TIF modelling requires the coupling of a dynamic, FD thermal reservoir model with FE fracture mechanics and a detailed wellbore system. TIF propagation changes are dynamically dependent on the following parameters [20]:

• Flow rate, pressure, and temperature of the injected fluid;
• Mechanical properties of the rock;
• Thermal properties;
• Reservoir properties; and
• In-situ stresses

The TIF modelling results reported in this study used a 2D FE fracture model coupled with a reservoir flow simulator and detailed wellbore model (see Figure 9). The shape of the TIF in the x-direction is shown in Figure 10.

The first step when modelling the development of a TIF requires estimation of the stress field in the reservoir. The subsequent changes in the stress field due to water injection are a function of the: (1) in-situ stress, (2) reservoir pressure, (3) reservoir temperature, (5) Poisson’s ratio, (6) Young’s modulus, and (7) thermal properties of the rock and fluid [20]. Once the stress field calculations are performed, the flow within the TIF is computed based on different TIF conductivity models. These conductivity models can either be dependent or independent on the TIF’s width. The leak-off flow can be determined from the difference between the pressures in the TIF and reservoir. TIF initiation and/or propagation is determined by a rock mechanical equation based on the stress–strain relationship. This relationship tests whether the TIF initiates and/or propagates based on the following parameters [21]:

• The difference between the BHP and minimum horizontal stress;
• Rock mechanical properties (e.g., Poisson’s ratio and Young’s modulus); and
• Critical stress intensity (K_IC) (i.e., rock toughness).

Based on the flow and rock mechanics equations, iterations are performed on the TIF geometry to calculate the stress intensity (K_I) at the tip of the TIF. The TIF will propagate only if the stress intensity calculated is higher than the input critical stress intensity factor (K_IC) (i.e., the fracture toughness of the rock). Figure 11 summarizes the different relationships used in the TIF model to examine TIF initiation and growth geometry within a dynamically coupled reservoir and wellbore system.

The TIF model used here employed has the following assumptions [22,23,24]:

• Formation rocks were continuous and linearly elastic;
• Mechanical properties of the rock were homogenous and isotropic;
• TIFs develop as planes (i.e., vertical TIFs were oriented orthogonally to the direction of the minimum horizontal stress);
• Reservoir porosity and permeability were independent of fluid pressure and formation stresses;
• TIF propagation was controlled by linear elastic fracture mechanics (LEFM) (i.e., K_IC was a measure of the intensity of the stress near the tip of the TIF required for TIF propagation to occur);
• Stress regime was normal faulting: σ_v > σ_Hmax > σ_hmin; and
• Magnitudes of the in-situ stresses varied with the depth.

The overall approach in this model was to subdivide the TIF into discrete elements (see Figure 10). The width, temperature, pressure and internal stress of the TIF were defined at the FE nodes. The detailed governing equations for the TIF model used here was comprised of the following [20,22]:

• Elasticity equations that related the pressure on the TIF’s tip to its width;
• Fluid flow equations that related the fluid flow and pressure in the TIF;
• Geomechanical equations that calculated the in-situ stresses, accounting for the effects of temperature and pressure changes on the reservoir stress field; and
• A TIF propagation criterion that related the stress intensity at the TIF tip to the critical stress intensity factor for the rock (K_IC).

The elasticity, fluid flow, geomechanical, and TIF equations are reviewed in the following sections.

4.2. Elasticity Equation

The elasticity equation solved the stress intensity (K_I) values at points around the TIF boundary. This required defined nodes with an identified pressure distribution within the TIF and stress over the TIF surface [20]. The TIF surface was defined with FE grid shape (see Figure 10), with the stress intensity values computed at discrete nodes on the TIF boundary. Stress intensity values at each node were obtained by computing the TIF’s width at all nodes of the XZ plane of the FE grid (see Figure 11). The TIF widths at different nodes and pressure were related by the following equation [20,21]:

$$\left(P_F-{\sigma }_{T,min}\right)\left(x,z\right)=\frac{G}{4\pi \left(1-\upsilon \right)}\int \left[\frac{\partial }{\partial x}\left(\frac{1}{R}\right)\frac{\partial w_f}{\partial x^{\prime }}+\frac{\partial }{\partial z}\left(\frac{1}{R}\right)\frac{\partial w_f}{\partial z^{\prime }}\right]d{x}^{\prime }{z}^{\prime }$$

(1)

where:

${\sigma }_{T,min}$	Total minimum horizontal stress normal to the TIF plane (y-direction in Figure 10) (psi)
$P_F$	Fluid pressure along the TIF
$w_f$	TIF width (ft)
$G$	Shear modulus (psi) defined by the Young’s modulus (E) and Poisson’s ratio ($\upsilon$). G is defined as:

$$G=\frac{E}{2\left(1+\upsilon \right)}$$

(2)

$\upsilon$	Poisson’s ratio
$E$	Young’s modulus (psi)
$R$	Distance between the source point (x’, z’) at which the integrand is evaluated and the field point at which the pressure is evaluated (x, z) (see Figure 10). R is defined as:

$$R=\left[{\left(x-x\prime \right)}^2+{\left(z-z\prime \right)}^2\right]$$

(3)

The details of the FE mesh and node generation as the TIF propagates, as well as the numerical method for finding an approximate solution for the TIF width (i.e., $w_f\left(x,z\right)$) can be found in [21,22].

4.3. Fluid Flow Equation

The fluid flow in a TIF is idealized as that of the laminar flow of an incompressible non-Newtonian fluid [21]. Here, the fluid flow in the TIF was assumed to flow between parallel porous walls. Leak-off through the TIF surface was determined by the difference between the pressures in the TIF and reservoir [22,24]. A 2D flow was obtained by integrating the governing equations through the width of the TIF. The flow and pressure within the TIF were related to the TIF’s leak-off rate by [20]:

$$-\int \frac{w_f{}^2}{12\mu }{\nabla }^2\left(P_F-\frac{\rho h}{144}\right)dV+\int M\left(P_F-P_p\right)dA+\frac{\int w_{f}dA-V_0}{\Delta t_f}-Q_f$$

(4)

where:

$P_p$	Pore pressure at far field (psi)
$\Delta t_f$	Time increment (sec)
M	Mobility connection factor, defined as:

$$M=8\frac{k_y}{\Delta Y}\sum \frac{K{r}_w}{{\mu }_w}$$

(5)

$k_y$	Reservoir permeability in the y-direction (see Figure 10)
$K{r}_w$	Water relative permeability (md)
${\mu }_w$	Water viscosity (cp)

The first term in Figure 5 is the flow rate within the TIF, the second term is the leak-off rate, the third term is the volumetric storage rate within the TIF, and the forth term was the total TIF injection rate added to the central node of the FE grid. The leak-off rate is assumed to be linear in the normal y direction (see Figure 10) from both vertical sides of the TIF. The mobility connection factor (M) in the leak-off term (see Equation (4)) was the mobility factor associated with the grid block intersecting the TIF with an area (A) and pressure difference between the TIF and pore pressure $\left(P_F-P_p\right)$ [20].

4.4. Geomechanical Equations

The stress calculation was an important part of coupling the fluid flow and elasticity equations within the FE grids of the TIF. The stresses that resulted from temperature and pressure changes in the reservoir were computed using the Goodier displacement method when a TIF was the absent (i.e., over the FD grids). Then, the computed stresses were interpolated into the 2D TIF surface using a simple model [20]. A single stress component normal to the TIF surface (i.e., the y-direction in Figure 10) was only considered when calculating the total stress, since it was necessary to calculate only the stress at the expected TIF surface.

4.5. Stress Calculation in the 3D FD Main Reservoir Grid

The stress field in the FD reservoir was calculated from the in-situ stress and rock mechanical properties. In-situ stress is defined as a function of the reservoir depth. In this study, the Goodier displacement method was used in the stress calculation [25]. The Goodier displacement method assumes zero displacement at the boundary of the reservoir model. This process was sufficient here, since vertical strain displacement and compaction were not considered in this research. The Thermo-elastic and poro-elastic stresses were calculated from the Goodier displacement potential using methods developed by Koning [24]. The stress resulting from temperature and pressure changes were described by the following equation [20,25]:

$${\sigma }_{T,min}={\sigma }_{hmin,i}+\Delta {\sigma }_T+\Delta {\sigma }_P={\sigma }_{hmin,i}+\Delta {\sigma }_y$$

(6)

where:

${\sigma }_{hmin,i}$	Initial minimum horizontal stress (psi)
$\Delta {\sigma }_T$	Thermo-elastic stress due to temperature change (psi)
$\Delta {\sigma }_P$	Poro-elastic stress due to pressure change (psi)
$\Delta {\sigma }_y$	Total stress reduction due to temperature and pressure change (i.e., sum of the thermo-elastic and poro-elastic stresses) (psi)

The total stress reduction due to temperature and pressure change was calculated by:

$$\Delta {\sigma }_y=\frac{\left(1+\nu \right)}{E}{\nabla }^2\mathrm{X}+(A_P\Delta P+A_T\Delta T)$$

(7)

$\Delta P$	Pressure difference between the flooded zone and reservoir pore pressure
$\Delta T$	Temperature difference between the injected fluid and reservoir (°F)
$A_P$	Poro-elastic constant (psi/psi), defined by:

$$A_P=\frac{\left(1-\frac{C_g}{C_b}\right)\left(1-2\nu \right)}{\left(1-\nu \right)}$$

(8)

$C_g$	Grain Compressibility (psi−1)
$P_F$	Bulk Compressibility (psi−1)
$A_T$	Thermo- elastic constant (psi/°F) defined by:

$$A_T={\alpha }_T\frac{E}{\left(1-\nu \right)}$$

(9)

${\alpha }_T$	Thermal expansion coefficient (1/°F)
${\nabla }^2\mathrm{X}$	Goodier displacement potential

The Goodier displacement potential was first computed over the entire FD grid [4]. The Goodier displacement potential (${\nabla }^2\mathrm{X}$) was given by the following equation:

$${\nabla }^2\mathrm{X}=\frac{\left(1+\nu \right)}{E}\left(A_P\Delta P+A_T\Delta T\right)$$

(10)

4.6. Stress Calculation for the 2D TIF Surface

The resultant change in the minimum horizontal stress ($\Delta {\sigma }_y$) due to temperature and pressure changes as computed using the Goodier displacement potentials were interpolated onto the 2D FE TIF surface once the TIF extended beyond its initial grid block [4]. The analytical model developed by Gonzales and Perkins [26] (see Equation (11)) with a circular shape function of 0.5 was used to the interpolate the stress from the 3D reservoir grids into the FE TIF surface [4]:

$$\Delta {\sigma }_y=0.5\left[A_P\left(T_i-T_R\right)+A_T\left(P_F-P_p\right)\right]$$

(11)

where:

$T_R$	Reservoir temperature (°F)
$T_i$	TIF surface temperature (°F)

4.7. TIF Propagation Criterion

Propagation of a TIF is controlled by the LEFM fracture criterion of. TIF propagation happens in such a way that the stress intensity factor (K_I) at each node is greater than the critical stress intensity factor (K_IC). Because the TIF width $(w_f)$ near the TIF tip region is proportional to the stress intensity factor at the boundary, the condition for TIF propagation can be given in terms of $(w_f$) as [4,22]:

$$w_f<w_c\mathrm{No}\mathrm{TIF}\mathrm{propagation}$$

(12)

$$w_f<w_c\mathrm{TIF}\mathrm{propagation}$$

(13)

where:

$$w_c=\frac{4K_{IC}\left(1-\upsilon \right)}{G}\sqrt{\frac{a}{2\pi }}$$

(14)

where:

$w_c$	Critical width at a fixed distance ($a$) from the tip of the TIF (ft)
$a$	Defined at a small distance from the tip of the TIF (ft1/2) (see Figure 4)
$K_{IC}$	Critical stress intensity factor ((rock fracture toughness)) (psi ft1/2)

The critical stress intensity factor ($K_{IC}$) is related to the additional pressure above the minimum horizontal stress required to open a TIF sufficiently to propagate. In this research, the critical stress intensity factor ($K_{IC}$) (i.e., fracture toughness) was obtained for brittle elastic solids from laboratory testing. Fracture toughness experiments performed on various rocks indicated that ($K_{IC}$) was of the order of 10³ psi-in^1/2 [22,27,28].

After computing $\left(w_f\right)$ from Equation (1), it was then compared to $\left(w_c\right)$ computed from Equation (14). Following the propagation criterion above, TIF propagation was then determined.

4.8. Solution Method

The FE elasticity and fluid flow equations within the TIF (i.e., Equations (1) and (4)) were combined and solved iteratively for the TIF widths $(w_f)$ using the Newton Raphson method. The iteration is performed in the following manner [8]:

• The combined FE equations solver was supplied by a total rate (Q_f).
• The geometry of the TIF was iterated until a constant TIF geometry was found (i.e., $w_f$ = $w_c$). The pressure at the center of the TIF (P_F) was returned and another level of iteration performed.
• The iteration continued until (P_F) and (Q_f) were consistent with the total well injection rate and BHP.

4.9. Implementation Workflow

The 2D FE TIF solution was coupled with the FD multiphase 3D flow in the reservoir, as shown in Figure 9. The reservoir flow and TIF modules were connected by a set of connection factors. The coupled elements included:

• Reservoir model: A multiphase 3D fluid flow (pressure and saturation) model created via an FD method and using black-oil pressure, volume, and temperature variables.
• Geomechanical solution: An FD method was employed to solve 3D stresses in the rock within the reservoir. The solution used the Goodier displacement potential method, which assumes zero displacement around the reservoir model (i.e., at the reservoir boundary). Once the stress was calculated for the 3D reservoir model, it was then interpolated into the 2D TIF model.
• TIF model: this was an FE grid with triangular internal and quadrilateral boundary elements (see Figure 11). The FE fluid flow and elasticity equations within the TIF were solved iteratively to obtain the TIF widths. The TIF propagation criterion was then applied to define the TIF shape and dimensions.

The coupling process and solution methods for all of these elements are shown in Figure 12 and Figure 13.

5. Sector Model Description

N Field Sector Model Description

The dimensions of the N Field sector model were 7874 ft × 6758.53 ft with a thickness of 903.8 ft. The 3D sector model grid system (see Figure 14) consisted of 26 × 22 × 474 cells giving an average size of 307 × 307.4 ft and a (fine) vertical thickness of 2 ft. The histograms for the porosity and permeability distributions of the cells are shown in Figure 15 and Figure 16.

The sector model had a total of four wells (see Figure 17). There were three horizontal wells (i.e., Producer NP4, Injector NI5_T, and Injector NI5_H) and one vertical well (Injector NI6). Producer NP4 and Injector NI6 were completed within the oil column, while Injectors NI5_T and NI5_H were completed within a (weak) aquifer partially masked by a shale layer. The reservoir also had a small gas cap. The vertical saturation map is shown in Figure 18.

6. Converting the Isothermal Reservoir Model to a Thermal Model

The isothermal model received provided fluid PVT properties at a range of pressures and single (reservoir) temperature of 224.6 °F (see Appendix A,

Table A1. Isothermal Fluid PVT Properties.

Temperature = 224.6 °F
Pressure (psi)	Bubble Point (psi)	Gas Oil Ratio (scf/STB)	Oil FVF RB/STB	Oil Viscosity cp	Gas FVF (RB/Mscf)	Gas Viscosity (cp)	Water FVF RB/STB	Water Viscosity (cp)
7251.89	738.415	1380.21	1.75093	0.283276	0.561812	0.033579	1.02306	0.32758
6890.03	738.415	1298.85	1.71203	0.292091	0.575612	0.032602	1.02428	0.32758
6528.16	738.415	1218.28	1.67351	0.303124	0.591175	0.031596	1.02551	0.32758
6166.31	738.415	1138.54	1.6354	0.316543	0.608871	0.030559	1.02673	0.32758
5804.45	738.415	1059.68	1.59769	0.332554	0.629183	0.029492	1.02795	0.32758
5442.58	738.415	981.727	1.56043	0.351413	0.652742	0.028393	1.02918	0.32758
5080.73	738.415	904.74	1.52362	0.37342	0.680391	0.027263	1.0304	0.32758
4718.88	738.415	828.775	1.4873	0.39893	0.71327	0.026103	1.03162	0.32758
4357.01	738.415	753.888	1.4515	0.428355	0.752957	0.024915	1.03284	0.32758
3995.15	738.415	680.158	1.41625	0.462165	0.801682	0.023706	1.03407	0.32758
3633.29	738.415	607.668	1.3816	0.500895	0.862681	0.022481	1.03529	0.32758
3271.43	738.415	536.518	1.34758	0.545137	0.940758	0.021254	1.03651	0.32758
2909.57	738.415	466.826	1.31427	0.59555	1.0433	0.020038	1.03774	0.32758
2547.71	738.415	398.739	1.28172	0.652833	1.18201	0.018856	1.03896	0.32758
2185.85	738.415	332.446	1.25002	0.717708	1.37632	0.017734	1.04018	0.32758
1823.99	738.415	268.179	1.2193	0.790879	1.66066	0.016701	1.0414	0.32758
1462.13	738.415	206.266	1.1897	0.872922	2.10233	0.015784	1.04263	0.32758
1100.27	738.415	147.181	1.16145	0.964077	2.85516	0.015003	1.04385	0.32758
738.415	738.415	91.6773	1.13492	1.06374	4.37215	0.014371	1.04507	0.32758
376.556	738.415	41.2191	1.1108	1.16895	8.84523	0.013896	1.04629	0.32758

). Extension of the fluid PVT properties to different temperatures was necessary to account for the temperature and pressure variations encountered during the productive life of the reservoir, due to injection of surface temperature water in the reservoir.

The methodology used was to match Table A1 to different correlations. The Glaso, Standing, Lasater, Vazques-Beggs, and Petrosky correlations for bubble point pressure, gas oil ratio (GOR), and oil formation volume factor (FVF) were all compared. The Beal et al, Beggs et al, Petrosky, and Bergman-Sutton correlations for oil viscosity were also compared. After checking the matching statistics for all correlations, the Vazques-Beggs and Beggs et al. correlations were selected.

The Vazques-Beggs correlations (see [29] and Table A1) contained equations for solving the GOR, oil FVF, and bubble point pressure. These were developed from data obtained from over 600 laboratory PVT analyses gathered from different fields around the world. These data covered a wide range of pressure, temperature, and oil properties [28]. The correlations used to convert the isothermal reservoir model into the thermal model are detailed in Appendix A.

7. History Matching the N Field Sector Model Without a TIF

After conversion into a thermal model, the N Field sector model was run and compared to historical data from Injector NI6. Injection Well NI6, being the point of interest of this study, was set to operate under THP control using the THP history data provided. Producer NP4 and Injector NI5_H controls were set to rate control using the historical data. The injection rate of Well NI6 was matched based on unreliable (estimated) reservoir pressure and BHP data. Any mismatch could be rectified later, after analyzing these uncertainties. The objective for running the simulation without TIF modelling was to answer the following questions:

• How did Injector NI6 perform during the entire period?
• When did the deviation from the historical data occur?
• Did the deviation correspond to the TIF occurrence period observed?

The resulting simulation and historical data are shown in Figure 19. The Injector NI6 simulation data showed a good match up to 520 days, after which increasingly greater deviations from the historical data occurred. This deviation corresponded to the TIF initiation expected to take place in the period between 520 and 700 days (i.e., between the dashed lines in Figure 19). The injection rate continued to increase during this TIF initiation period, even though the BHP was still decreasing. This was especially clear after the 700 days points (i.e., the second dashed line in Figure 20). This was a clear indication of TIF occurrence. Modelling the injection rate with a radial outflow equation did not produce a match to the reliable data after 520 days. Thus, it could be confidently concluded that TIF initiation occurred after 520 days.

8. In-Situ Stress

The in-situ stresses values for the N Field and rock mechanical properties had to be matched because these are fundamental parameters governing TIF occurrence. TIF initiation and growth was observed in Well NI6 above, and by analogy probably also occurred in other N Field injectors. Well tests and well logs were available for Well N-8. The location of Well N-8 relative to Injector NI6 is shown in Figure 21

• The direction of the principle in-situ stress in the “N” field is vertical
• The stress regime was normal faulting: σ_v > σ_Hmax > σ_hmin.
• All measurements and estimations of in-situ stress and rock mechanical properties for Well N-8 were assumed to also apply to Injector NI6.

8.1. Estimation of Vertical Stress

The magnitude of the vertical stress (σ_v) equal to the weight of the overburden was calculated from the integrated bulk density log for Well N-8. Figure 22 shows the bulk density (ρ_b) from the top of the logged interval to the top of the formation. The weight of the overburden for N-8 was calculated down to a depth of 2970 m (9744 ft) (i.e., the perforated interval in Well NI6). There were different rocks below the seabed. The average densities of these rocks at various depths can be found in Figure 22. These values, in SI units, were eventually converted to oilfield units, 1000 kg/m³ for 130 m (i.e., seawater: 1850 kg/m³ up to 1000 m, 2400 kg/m³ up to 2300 m, and 2700 kg/m³ up to 2970 m). These values were then used to estimate the vertical stress (σ_v):

$${\sigma }_v={{\displaystyle \int }}_0^z\rho \left(z\right)gdz$$

(15)

$${\sigma }_v=\mathrm{Weight}\mathrm{of}\mathrm{Sea}\mathrm{Water}+\mathrm{Weight}\mathrm{of}\mathrm{Rocks}$$

(16)

• ${\sigma }_v=9.83\left[\left(1\times 130\right)+\left(1850\times 870\right)+\left(2400\times 1300\right)+\left(2700\times 670\right)\right]$
• ${\sigma }_v=6.56\times {10}^{7}Pa=9514\mathrm{psi}$

A vertical stress value of 9514 psi was calculated, transforming it into a stress gradient of 0.97 psi/ft. This vertical stress gradient was used to calculate the vertical stress during TIF modelling in the remainder of the analysis.

8.2. Magnitude of the Minimum Horizontal Stress

One of the most common methods in the literature for determining the magnitude of the minimum horizontal stress is to use leak off tests (LOTs). Oyeneyin [30] indicated that the LOT method was not the most accurate because such tests are not performed in reservoir intervals and can only indicate a stress between the vertical and minimum horizontal stresses. However, the best measurement available was via LOTs performed on Well N-8 at the two different intervals used in this study (i.e., 824 and 2678 m).

The LOT at 824 m true vertical depth (TVD) in Figure 23 showed an equivalent mud weight of 13.5 ppg (or 0.7 psi/ft) at the leak-off pressure. A higher LOT value of 14.56 ppg (or 0.76 psi/ft) was recorded at 2678 m TVD (see Figure 24). The average of these two values (0.73 psi/ft) was used as the starting point in the history matching analysis.

8.3. Magnitude of the Maximum Horizontal Stress

The magnitude of the maximum horizontal stress (σH_max) was the most difficult stress to estimate. However, drilling experience and wellbore imaging offer important information regarding the magnitude and orientation of (σH_max) [31]. An approximate value for (σH_max) can be obtained using a stress polygon if the stress regime, vertical stress from the density log and minimum stress from the LOTs are all known.

A stress polygon defines all possible maximum and minimum stress magnitudes according to Anderson’s theory and Coulomb faulting theory [31]. The limiting ratio of maximum and minimum principle effective stress is given for normal faulting by:

$$\frac{{\sigma }_v-P_P}{{\sigma }_{hmin}-P_P}\le [{\left({\mu }^2+1\right)}^{\frac{1}{2}}+\mu ){]}^2$$

(17)

Here, $\mu$, the coefficient of friction, was assumed to be 0.6, and Pp, the pore pressure, was assumed to be a gradient of 0.433 psi/ft. Figure 25 is the resulting stress polygon at a depth of 8,786 ft, the LOT depth.

In principle, the Figure 25 stress polygon permits a wide range of (σH_max) values at a depth of 8786 ft TVD. The regional stress regime was normal faulting, which when combined with the deep LOT gave a narrower range of possible stress values. The average value (7513 psi or 0.86 psi/ft) was chosen as the starting point for the history matching analysis.

8.4. Orientation of the Minimum Horizontal Stress

A TIF propagates perpendicular to the direction of the minimum horizontal stress (σh_min). The formation micro imager (FMI) log often provides information regarding stress direction. Well N-8’s FMI did not show any failure features that could have been caused by drilling. However, Well N-20’s FMI (see Figure 26) clearly showed tensile fractures propagating in an E–W direction, hence, (σH_min) had a N–S orientation.

9. Rock Mechanical Properties

An estimate of the rock mechanical properties for the reservoir were required for TIF modelling. The sonic log for Well N-8 was used to derive the Young’s modulus and Poisson’s ratio, since no measurements were made of the core samples. Fracture toughness and Biot’s coefficient were estimated from published correlations.

9.1. Young’s Modulus and Poisson’s Ratio

The dynamic Young’s modulus (E) and Poisson’s ratio (v) were derived from the log-based compressional ($V_p$) and shear wave $\left(V_s\right)$ velocities [32]:

$$v=\frac{\frac{1}{2}{\left(\raisebox{1ex}{$V_p$}\!\left/ \!\raisebox{-1ex}{$V_S$}\right.\right)}^2-1}{{\left(\raisebox{1ex}{$V_p$}\!\left/ \!\raisebox{-1ex}{$V_S$}\right.\right)}^2-1}$$

(18)

$$E=\rho \times V_S{}^2[\frac{3V_P{}^2-4V_S{}^2}{V_P{}^2-4V_S{}^2}$$

(19)

Compressional $\left(V_p\right)$ and shear wave $\left(V_s\right)$ velocities were estimated using Well N-8 log measurements (see Figure 27). The compressional $\left(V_p\right)$ and shear wave $\left(V_s\right)$ velocities were estimated to be 4545 and 2702 m/s, respectively. The Young’s modulus (E) was determined to be 47.04 GMpa (6.82 × 10⁶ psi) and the Poisson’s ratio (v) 0.23. These values were used as a starting point in the initial history matching analysis.

9.2. Fracture Toughness

The TIF propagation criterion employed in this research was the fracture toughness (K_IC) as a measure of the rock’s resistance to fracturing. A fracture will propagate once the stress intensity factor of a Type 1 fracture (K_I) is sufficiently larger than (K_IC) [33]:

$$K_I\ge K_{IC}$$

(20)

(K_IC) can be measured by a laboratory experiment or field scale fracture test. The statistical relationship between (K_IC) and the dynamic Young’s modulus (E) for sandstone developed by Zhix [33] was used becasue core measurement data were not available:

$$K_{IC}=0.0215\times E+0.2468$$

(21)

Figure 21 gave a fracture toughness of 1.26 Mpa/$\sqrt{m}$ (330 psi/$\sqrt{ft}$ ). This value was used as a starting point in the history matching analysis.

9.3. Biot’s Coefficient

Biot’s coefficient was used to estimate the reduction in total stress due to reservoir cooling by the cold injection water. The Wu empirical correlation for consolidated sediments was used here [32]:

$$\alpha =1-{\left(1-\varnothing \right)}^{3.8}$$

(22)

where:

• $\alpha$ Biot’s coefficient
• $\varnothing$ Porosity

The perforated formation’s porosity (see Figure 27) was 0.18 at a TVD of 9744 ft, giving an estimated Biot’s coefficient of 0.53.

10. Uncertainty Analysis

The initial estimates of the NI6 well data discussed above were used for the initial TIF performance model. It was understood that these estimates would be accompanied by a great deal of uncertainty. This section assesses the uncertainty related to the NI6 reservoir pressure. TIF growth is also evaluated. The in-situ stresses and rock mechanical data were input into the reservoir simulator coupled with the geomechanical solution and TIF model. Injector Well NI6 was operated under THP control using the historical data, while Producer NP4 and Injector NI5 were rate controlled to reproduce the field history data. The average surface temperature of the NI6 injection water during the historical period was 70 °F. The simulated “with TIF” and the modelled “no TIF” results were in good agreement, (see Figure 28) up to 520 days. Section 7 showed that TIF initiation occurred between 570 and 700 days

The with TIF simulation predicted that TIF initiation occurred at 620 days. The TIF initially delivered an injection rate that exceeded the historical data, but it rapidly decreased to a much lower no TIF rate after ~750 days. This divergence between the simulated and historical data might be due to the use of incorrect values for reservoir pressure, injected water temperature or thermal and geomechanical rock properties. This notion is evaluated in the following sections.

10.1. Reservoir Pressure

Average reservoir pressure data for Well NI6 was measured by long-term fall off-tests during the 2009–2012 period (see

Table 1. Average Reservoir Pressure Measured by Long-term Fall off-Tests for Well NI6.

Well NI6 Shutdowns		Average Regional Pressure (psi)
After 266 days	July -09	4480.5
After 702 days	Oct -10	4524
After 1277 days	May -12	4263

). Figure 29 shows that the TIF simulated reservoir pressure around Well NI6 did not match these fall-off test measurements. The simulation initially overestimated the actual reservoir pressure by 100 psi, increasing to >400 psi after four years. The simulated reservoir pressure needed to be matched to reproduce the fall off-tests that measured the actual reservoir pressure.

10.2. Wellbore Modelling

Figure 30 shows that the NI6 wellbore model did not match the historical data when the well was controlled by the wellhead pressure. Section 11.2 describes the calibration of the wellbore model to reproduce the BHP history.

10.3. Type and Surface Temperature of the Injected Water

The injected water’s temperature at the surface as well as its type (or properties), are important aspects of TIF modelling. The historical data for the surface water temperature for Injection Well NI6 (see Figure 31) showed variations between 60 and 90 °F that were difficult to properly include in the reservoir simulator. Knowledge of the type of water injected and its properties was required, in addition to its temperature. Clean (i.e., low solids) sea water (SW) typically has a lower temperature and creates a higher conductivity TIF, while produced water (PW) (i.e., with a higher solid content) leads to TIFs with lower conductivity and higher propagation pressure.

11. History Matching with the TIF Modelling Workflow

The objective here was to match the history of a key well (i.e., NI6) that showed a clear indication of a TIF. History matching with TIF modelling can be a time-consuming process for two primary reasons. First, history matching is performed manually due to the dynamic nature of TIFs and persistent transmissibility changes that occur with time. Second, history matching involves geomechanical aspects (i.e., in-situ stresses and rock mechanical properties) in addition to reservoir properties and wellbore modelling. This complex history matching problem was divided into a two-step process:

• Figure 32 indicating the workflow used to achieve a good match of the no-TIF simulation and historical data. The injection rate of Well NI6 was not matched at this stage. This workflow is discussed further in Section 11.1 and Section 11.2.
• Figure 33 shows the workflow used to obtain a good history match after a TIF has occurred. The in-situ stresses and rock mechanical properties were matched here. A good match to the injection rate of NI6 was the ultimate goal at this stage. This workflow is discussed further in Section 12.

11.1. Global History Matching

There were multiple parameters to be matched. A top-down approach from the field level to the well level was adopted to minimize the effort required to obtain the final history match.

Regional Reservoir Pressure

Table 1 provides three data points for matching Well NI6’s regional reservoir pressure. Reservoir pore volume and compressibility were modified at this stage to match the simulated reservoir pressure to the Table 1 pressure fall-off data. The final matched reservoir pressure is shown in Figure 34.

11.2. Local History Matching

Well NI6 BHP

A detailed wellbore model was built in PROSPER^TM i.e., well performance and design software developed by Petroleum Experts based on the well’s completion schematic (see Figure 35), a survey, BHP measurements during an injection test, etc. These data were used to construct a mathematical model (i.e., wellbore model) in PROSPER that was calibrated with the measured data. The following workflow was applied:

• A critical review of the raw well test data (see Figure 36) addressed how reliable the measurements were and how the test data compared to the historical trends.
• Three test points from the injection test (see Figure 36) were selected for further analysis
• All vertical lift performance (VLP) correlations in PROPSER were compared and the Petroleum Experts 2 correlation was selected since it reproduced each well test with reasonable accuracy (see

  Table 2. Difference between Measured Test Data and Calculated Using PROSPER.

  Test Points	Time (min)	Liquid Rate (STB/d)
  		Measured Values	Calculated Values	% Difference
  Point 1	60.3	10,860	10,409	−4.33
  Point 2	133.6	5013	5045	0.63
  Point 3	178.2	19,215	17,968	−6.94

  ).
• Finally, the simulated BHP successfully reproduced the historical BHP (see Figure 37).

12. Final Matching with TIF Modelling

The reservoir pressures and BHPs related to Well NI6 were matched and the focus directed to the thermal and geomechanical properties. The injection rate of Well NI6 was then matched with the well history by following the workflow outlined in Figure 32.

12.1. Injected Water Temperature

Before going into detail regarding the modified geomechanical and thermal properties, it is important to discuss other uncertainties (e.g., water temperature and type). The operator confirmed that only SW was the source of the injected water. PW had not been reinjected. The field was located near deep water and distant from the coast. The SW could thus be assumed to have a low solids content and therefore a high TIF conductivity that would not affect TIF propagation. The field data did not show a clear correlation among the injection temperature, injection rate, and TIF after 690 days (see Figure 38). The injection of colder water has been found to improve a well’s injectivity index in several published field histories [1,3,34]. Further theory supports the notion that cooling the formation affects the TIF dynamics by encouraging TIF growth and greater TIF conductivity. However, an average water injection temperature of 70 °F had to be assumed, due to the technical limitations of the geomechanical simulator. The TIF was assumed to have an infinite conductivity based on the above assumption that water quality was not an issue for the N Field.

12.2. Geomechanical and Thermal Properties

History matching of the injection rate without TIF after correcting the reservoir pressure for Well NI6 is shown in Figure 39. It should be noted that the simulation deviated from the historical data after 690 days. The In-situ stresses, rock mechanical properties, and thermal properties were then modified from values determined in Section 8 and Section 9 to match the rate injection history using the workflow outlined in Figure 32. The final history matching with TIF modelling is shown in Figure 40.

Matched geomechanical and thermal properties are shown in

Table 3. Initial and Final Matched Geomechanical and Thermal Properties.

Property Type	Property Name	Initial	Matched
Rock Mechanical Properties	Young’s modulus (MMpsi)	6.8	2
	Poisson’s ratio	0.23	0.35
	Fracture toughness (psi/ft1/2)	330.48	150
	Biot’s coefficient	0.53	0.67
In Situ Stresses	σv (psi/ft)	0.97	0.97
	σHmax (psi/ft)	0.86	0.86
	σhmin (psi/ft)	0.73	0.65
Thermal Property	Rock thermal expansion coefficient (1/F0)	3.20 × 10−4	1.20 × 10−5

. Based on the Table 1, the most uncertain parameters were the Young’s modulus, Poisson’s ratio, fracture toughness, and the minimum horizontal stress. The matched Poisson’s ratio was higher than the original, estimated value (see Section 9.1), whereas the final Young’s modulus was lower than the calculated value. These two elastic coefficients were used to calculate the stress reduction (or thermoelastic stress) due to cold water injection. Therefore, these values were used to calculate the reduction in minimum horizontal stress and, hence, influence the initiation and propagation of the TIF. The initial minimum horizontal stress value was reduced during the history matching process (it had been overestimated in Section 8.2). It was necessary to reduce this stress in order to obtain a good final historical match. Fracture toughness was also reduced to a lower value, since the initial estimate did not allow for fracture initiation or/and propagation. The model was finally used with confidence to test alternative development scenarios for injection well completion.

There were many sources of uncertainty in this study. Specific laboratory tests on core data were not available, hence; generalized correlations were employed. These values were also assumed to be isotopic and homogenous, which is far from the reality in complex reservoirs. The geomechanical simulator was limited to a constant average value for the injected water temperature for the entire injection period.

12.3. TIF Growth Dynamics

TIFs were initiated and propagated based on the modelling analysis of the first 690 days and onward. This time corresponded to the sudden sharp decrease observed in the Well NI6 BHP history, as shown in Figure 37. It was also noted from the simulation (see Figure 41) that the TIF geometry propagated on both sides of NI6 and parallel to Producer NP4, as expected from TIF direction analysis. The TIF half-length increased sharply within a few days, as shown at Point 1 in Figure 42, and then stabilized for some time. The TIF then propagated and increased after more cooling took place near the TIF tip at Point 2 in Figure 42. The TIF half-length then slightly decreased at point 3 in Figure 42. This decrease was not supported by the physical changes observed and could be due to numerical error.

Even though the TIF half-length could not be correlated to the injected water temperature (since it was assumed as one value 70 °F), the cooling effect can clearly be identified in Figure 42. The TIF tip stopped propagating when the pressure inside the TIF was not high enough to propagate the TIF (Point 1). After the rock cooled further over time, the TIF propagated until it reached warmer rocks (Point 2). Even though the TIF half-length was not certain and could not be confirmed with data, it physically behaved as expected.

13. Conclusions

History matching of a dynamic phenomenon like TIF formation can be a difficult task, due to the interaction of various complex processes. The workflows developed in this research and the associated coupling approach are new and can help to address different challenges during waterflooding operations.

In this study, the performance of Injector NI6 was modelled without a TIF and history matched to the injection data. The magnitudes and orientations of the in-situ stresses were determined, along with the rock mechanical properties, based on the data provided, well logs, and published correlations. These geo-mechanical properties, together with upgrading to a thermal PVT model, allowed for modelling and history matching of Injector NI6’s performance with a TIF. It was concluded that:

• The direction of maximum horizontal stress, and TIF propagation was most likely E–W, with a normal faulting regime.
• The surface water temperature and both the volume of the injected water and TIF propagation, were well-correlated, especially after the onset of the TIF at 690 days.
• History matching with TIF modelling proved to be time-consuming, due to the dynamic nature of TIFs and the presence of geomechanical uncertainties.
• Two workflows were developed for history matching for Well NI6 data. The reservoir and well parameters were first matched to the field data prior to TIF initiation at 690 days. The geomechanical properties were adjusted by history matching the subsequent field data to the TIF model.
• The model results without TIF were matched to the NI6 injection rate data before day 690. Inclusion of the TIF and suitable modifications to the initially estimated values of the geomechanical properties were required to obtain a comparable match after this date.
• The reservoir and geomechanical models were improved and validated once the modelled and observed TIF onset and propagation periods were confirmed
• The resulting history-matched N Field sector model can be used with confidence in future studies.