# A synthetic example

In this chapter, I formulate a new equation to invert for pressure and saturation changes in a thick, compacting, chalk reservoir. This is achieved via synthetic modelling of changes in different dynamic properties such as pressure, saturation and compaction, using the fluid flow simulator ECLIPSE and then translating them into elastic properties using a petro-elastic model. Subsequently, by employing backward engineering, I decompose the composite relative impedance to analyse its individual components and workings in detail in order to then recreate the composite relative impedance using a simpler approximation or a proxy model. Physical phenomena such as water weakening and compaction notorious in chalk reservoirs are accounted for in the proxy model, and also described in this chapter. The validation of the proxy model is carried out in both forward modelling and inversion, and the benefits of incorporating constraints such as engineering concepts and bounding values are also discussed at length.

Chapter 1 introduced various quantitative approaches in estimating the reservoir dynamic properties such as pressure and saturation changes, specifically the use of data-driven proxy models (MacBeth et al., 2006, Floricich et al., 2006, Corzo et al., 2013, Alvarez and MacBeth, 2013). These models relate the changes of seismic amplitude to the changes in dynamic properties. The inspiration to break away from these map-based models arises from their complicated application on the Ekofisk field, which is a thick, multi-cycle, compacting chalk reservoir. This chapter describes the derivation of a more appropriate proxy model, through acoustic impedance and time-shifts.

## 1.1 Simulation Model Description

We can study fluid movements, pressure propagation rates and drainage patterns using a numerical reservoir model and simulations to evaluate changes in a reservoir during production. I will first describe in detail the characteristics of the Ekofisk reservoir simulation model which will be used to generate synthetic elastic properties to help us to better understand how pressure, saturation and compaction affect the elastic parameters. This model is unique in that it represents the complex dry compaction and water weakening mechanisms that greatly influence how reservoir fluids are produced, injected and affect the rock frame.

The simulation model has 128 x 155 x 22 (436480) grid cells, with an average cell size of 100 x 110 x 45 m corresponding to an area of approximately 11 km x 6 km. The model is built to be consistent with the geological features, flow units and fault planes; hence, the grid uses non-vertical pillars and irregular cells (corner-point geometry (CPG)). The reservoir model has an anticlinal structure. The reservoir is composed of chalk material; its initial porosity ranges between values of 0.02 to 0.48 and its horizontal permeability ranges between 0.0007 to 2000mD and in the vertical direction ranges between 0.00007 to 200mD. The static parameters are well constrained by the large amount of well data. Data from wells and an acoustic impedance background trend were used in the porosity modelling. There are 23 geological horizons and 22 geological layers in the model. The preferential fluid flow patterns and, therefore, the preferential changes in seismic attributes are influenced by the rock types of the model. The simulation model with all 22 layers of initial porosity is illustrated in Figure 4.1. The fluid properties and initial reservoir conditions are shown in Table 4.1.

The Ekofisk field has a large variability in reservoir quality. The fracture network indicators are characterised deterministically by large amount of data using both static (e.g. fracture distribution) and dynamic observations (e.g. well tests, water breakthrough). The effective permeability model is a combination of two properties: fracture-enhanced matrix permeability and fracture network k permeability based on 14 different indicators: interpreted fractures, flow path analysis, distinct water breakthrough, tracer, super tracer, temperature anomaly, pressure supported data, interference test, rapid gas-oil ratio (GOR) increase, mud loss, PLT kick, fracture, fracture cluster and observations from 4D seismic data. These fracture network indicators are mapped directly into the flow simulation grid. This deterministic 3D mapping forms the basis of a high-contrast single porosity and permeability model, which is used as input for dynamic simulation (Tolstukhin et al., 2012).

##### Figure 4.1: The fluid flow simulation for the Ekofisk field with 128x155x22 grid cells in total. High porosity regions are found in the crest of the reservoir.

Property |
Description |

System | Gas, oil, water, dissolved gas in live oil, vaporized oil in wet gas |

Oil gravity (API) | 38 |

Initial reservoir pressure (psia) | 7150 |

Temperature (^{o}F) |
268 |

Gas gravity | 0.73 |

Salinity (ppm) | 65,000 |

Gas density (lb/ft^{3}) |
0.055 |

Oil density (lb/ft^{3}) |
52.23 |

Water density (lb/ft^{3}) |
62.37 |

Residual oil saturation (S_{orw}) |
0.275 |

Connate water (S_{wi}) |
0.05 |

Residual gas saturation (S_{org}) |
0.04 |

Critical gas saturation (S_{gc}) |
0.36 |

#### Table 4.1 Fluid reservoir properties and initial reservoir conditions used in the reservoir simulation.

### 1.1.1 Rock typing

The modelled porosity, effective permeability and fracture index are used to calculate the eight different rock types in the simulation model, where the function of rock typing is to allocate relative permeability and saturation region information for each of the cells in the model. The model is a single porosity model; hence the fracture or matrix exchange needs to be covered by pseudo relative permeability functions (ConocoPhillips Internal Report). The absolute permeability is a property of the reservoir porosity and is a measure of the capacity of the rock to transmit fluids. When two or more fluids flow at the same time, the relative permeability of each phase at a specific saturation is the ratio of the effective permeability of the phase to the absolute permeability (Ahmed, 2010). The same rock type can have varying porosity and permeability but correspond to the same relative permeability curves. The mapping criteria for the eight different rock types are shown in Table 4.2 below:

Rock Type |
Effective permeability |
Porosity Description |
||

20 |
>100mD | Thief zone | No collapse | |

18 |
Collapse | |||

16 |
>15mD | Fractures | >32% | High porosity |

15 |
<32% | Low porosity | ||

14 |
5mD> & <15mD | Intermediate | >32% | High porosity |

13 |
<32% | Low porosity | ||

12 |
<5mD | Matrix | >32% | High porosity |

11 |
<32% | Low porosity |

#### Table 4.2 Rock types in the Ekofisk field are divided based on the effective permeability and porosity.

Major fracture corridors along the faults with effective permeability greater than 15mD are categorised under rock types 15 and 16, whereas the background matrix is mostly described by rock types 11 and 12. Rock types 18 and 20 are thief zones, which are horizontal permeability conduits which contribute to early water breakthroughs in producer wells. The rock types are shown in Figure 4.2, where majority of the cells are described as matrix.

##### Figure 4.2: A map view of the Ekofisk simulation model (layer 11) coloured by rock types. Fractures are highlighted as cyan.

### 1.1.2 Geomechanical changes handled by reservoir simulation

The water weakening phenomenon was widely studied in the 90s, with many notable publications (Newman, 1983, Schroder and Shao, 1996, Delage et al., 1996, Risnes and Flaageng, et al., 2004, Risnes at al., 2005, Austad et al., 2008). The change in pore volume has been shown in laboratory experiments where the rock compacts as a function of increase in effective stress, rocks that are more porous also underwent more dramatic compaction gradient than less porous chalk. According to the compaction model (Sylte et al., 1999) provided by the operator (ConocoPhillips), porosity rebound does not occur due to unloading events. This means the rock behaves in an inelastic manner, even after pore pressure has increased to its initial condition. The same relations also apply to the relationship between permeability and stress: for less porous chalk, there is no change in fracture conductivity or matrix permeability as effective stress reduces in the reservoir. Porous chalk will undergo a more severe reduction in permeability in the presence of water coupled with pressure depletion compared to only pressure depletion. The dry compaction and water weakening mechanisms are incorporated into the simulation by using the keyword ‘ROCKCOMP’ in ECLIPSE. The compaction tables are read into the simulator using the ‘ROCKTAB’ keyword.

## 1.2 Dry Compaction and Water Weakening

Dry compaction is where the rock reduces in pore volume when the fluid pressure falls or the effective stress increases. This behaviour can have both positive and adverse effects on the reservoir, such as adding significant energy to the reservoir but also causing massive compaction which translates into subsidence of the seafloor. However, some rocks, typically chalks, will exhibit additional compaction when the water meets oil bearing rock, even at constant stress. Laboratory work by Newman (1983) and Loe et al. (1992) on North Sea chalk showed immediate and dramatic weakening of water-free chalk when injected with sea water. From the work of Smith et al. (2002), based on an extensive database of chalk core measurements, it has been shown that the pressure drawdown mechanism or dry compaction creates a gentler compaction trend compared to water weakening, at the same effective stress. It also showed that both compaction mechanisms are weaker for chalk with lower initial porosity.

Laboratory tests also indicate that the relative amount of water weakened chalk is a linear function from zero water saturation to the critical water saturation, which is when the entire matrix is fully water wet and the chalk is fully water weakened. The critical water saturation or the maximum attainable water saturation (

Swmax

) has been determined as 0.325 in the Ekofisk field (Smith et al., 2002), although it was initially reported to be 0.25 (Sylte et al., 1999). The rock is known to be ‘fully-wetted’ when the maximum water saturation is reached, as shown in Figure 4.3. The chalk does not need to be saturated at 1-

Sorw

(residual oil saturation) to be consider fully wetted. The fully wetted stage is when the chalk undergoes maximum compaction after water invasion from 0 to 33% at a given effective stress for a specific porosity. The maximum attainable water saturation (

Swmax

) is unique to the water weakening phenomenon. The deformation mechanisms described here are paramount in capturing the changes of porosity as a function of effective stress and changes in fluid composition. The effects of water weakening are first visible in samples of oil saturated chalk appear to be 2-3 times as strong as water saturated samples (Risnes and Flaageng, 1999), as illustrated in Figure 4.4. This further complicates the stress-strain behaviour of chalk, as samples with the same porosity, will display different curves for different water saturations (shown in Figure 4.5).

## 1.3 Separating Pressure and Saturation changes in a thick versus thin reservoir

In this section, I outline how the thickness and heterogeneity of the reservoir determine techniques carried out in 4D seismic data interpretation, and more importantly, in efforts to separate pressure and saturation changes. I will also explain some of the advantages of working in the impedance domain instead of the reflectivity domain. Figure 4.6 shows a spectrum of the different types of reservoir from ultra-thin to thick, heterogeneous, compacting reservoirs. My field of interest, Ekofisk is categorised under thick, heterogeneous and compacting reservoirs, similar to reservoirs such as Luconia, Sleipner and the Dan field.

##### Figure 4.6: Different reservoirs categorised based on reservoir thickness and heterogeneity.

Two synthetic examples are generated to demonstrate why elastic properties are the most suitable attributes to decompose pressure and saturation signals for a thick and heterogeneous reservoir. In the first case, shown in Figure 4.7, I modelled a thin reservoir with an injector injecting water into an oil leg, and observed an intermixture of pressure and saturation changes in the compartment. The change in pressure is assumed to be 1000psi and the change in water saturation is 30%. The reservoir is thin; hence the gravitational effects of fluids are less apparent. As shown in Figure 4.7 (a) and (b), pressure diffuses across the entire reservoir and water slumps due to the gravitational effect. Assuming there is a producer up dip of the injector (a typical scenario), the pressure gradient is shown in Figure 4.7 (b), and water is shown to have migrated towards the producer, potentially with water coning around the producer. The pressure effect is more dominant than saturation changes; therefore, an overall decrease in impedance of -3% is observed at the injector based on rock physics modelling.

The seismic profile in Figure 4.7 (c) at the injector well shows the quadrature phase (-90^{o}) amplitude of the baseline (the blue trace) and of the monitor (the red trace). The dashed black line is the difference in amplitude of the quadrature phase before warping, and the solid black line is the time aligned quadrature phase difference in amplitude, which shows a negative value. In the case of any zero-phase data, the energy peaks at interfaces, which are at the top and bottom of the reservoir. The reason I convert seismic difference to quadrature phase (-90^{o}) is to display the energy difference within the reservoir interval, and it is also useful for facilitating volume-based interpretation techniques (Johnston, 2013).

From a map view of a seismic attribute such as the difference in root mean squared amplitude shown in Figure 4.7 (d), it can be seen that both top and base of the reservoir show a softening response and it will not be possible to separate these changes easily. If the reservoir is thinner than tuning thickness, the base of the reservoir might not be resolved. If there is prior knowledge that the injection is in an oil leg and not the water leg, such that the depth of the original oil-water contact is known, it can be inferred that the softening is attributed to a combined effect of both pressure and saturation. Moreover, if the reservoir corresponds to a half cycle of the seismic amplitude, these softening responses can be mapped spatially. In the case in point, I cannot proceed further in terms of separating pressure and saturation solely by analysing the amplitude difference, unless I invert the amplitude difference to impedances and subsequently carry out a rock physics transform to separate the individual dynamic properties.

With a thick and heterogeneous reservoir, for which there is potentially a vast amount of information, it is apparently more complicated to interpret the seismic anomalies, and the same anomalies can be explained by many different production mechanisms. In Figure 4.8 (a), using the same example as the thin reservoir with an injection into the oil leg, a different anomaly is observed. I assumed a greater gravitational effect in the thick reservoir, where a longer transition zone exists between oil and water, as shown in Figure 4.8 (b). Again, assuming a producer up dip of the injector, we expect to see slumping of the heavier reservoir fluid (water) at the injector well, and water coning at the producer well. Of course, the coning is subjected to production rate at the well. On the other hand, we expect pressure to diffuse uniformly throughout the reservoir interval. The top reservoir therefore will have a stronger pressure signal compared to the bottom of the reservoir, because the water coning effect is stronger at the base of the reservoir. In the thick, heterogeneous reservoir, I can observe both the changes in pressure and saturation, spatially, and also the vertical extent of the injection. This is due to the varying impedance contrast of the different chalk layers, which is an intrinsic property of the rock, like porosity and permeability. Figure 4.8 (c) shows the quadrature difference as a solid black line. A softening signal (negative amplitude difference) near the top of the reservoir is observed, whereas a hardening (positive amplitude difference) is observed at the base of the reservoir. Both of these signals are highlighted by the red and blue solid arrows. These are genuine signals caused by dynamic changes. The hardening and softening responses in the middle of the reservoir are due to side-lobe interference.

In multiple stacked reservoirs similarly to those demonstrated in this synthetic example, the side lobes can interfere with or be confused with the primary difference signal. The solution for such quadrature shortcomings is 4D seismic inversion. Although I cannot accurately separate the pressure and saturation signals within the reservoir completely, it helps us to interpret the vertical extent of the water propagation, based on the polarity of the amplitude difference. More sophisticated methods, such as impedance inversion will be required to ultimately separate both of these effects. If I tackle this problem from a map-based method, as illustrated in Figure 4.8 (d), different responses will be seen at distinct interfaces of the reservoir. The separation of pressure and saturation is case-dependent with this type of reservoir, such that the top reservoir and perhaps the intra-reservoir layer will yield softening responses due to pressure and the base reservoir shows hardening due to water saturation. This interpretation strategy is also applied in the Andrew field. From the 4D difference cross-section in Figure 4.9(a), the extent of gas and water changes can be visually inferred from the polarity alone. 4D water and gas migration maps were computed by summing the positive and negative differences across the oil column (Trythall et al., 2003). The map generated for the top reservoir in Figure 4.9(c) shows gas saturation changes based on the softening signals and vice versa for the map generated from the base of the reservoir, shown in Figure 4.9(b), where hardening is correlated to water injected.

In the case of Sleipner, due to the nature of the reservoir with many intra-reservoir shales in the Utsira sandstone, the extensive propagation of gas was effectively revealed by the shales when the gas was trapped beneath thin shale layers. Thus, the thick, intra-shale reservoir enabled the successful qualitative monitoring of the migration of carbon dioxide gases (Figure 4.10). Thus, there are also many advantages to a thick, heterogeneous reservoir and careful interpretation can help us to gain insights into the dynamic changes of the reservoir without impedance inversion, to a qualitative extent. However, due to the interference of side-lobes and without prior knowledge regarding the geology, it will be difficult to interpret the seismic anomalies, it is henceforth beneficial to use a layer property such as impedance instead of an interface property like amplitude.

## 1.4 Pressure and Saturation Sensitivity on Elastic Properties

In this section, I will present the responses in various elastic properties influenced by primary production mechanisms in the Ekofisk field, modelled for the LoFS surveys. This sensitivity analysis was carried out to understand the different production mechanisms in isolation, particularly for pressure, saturation and porosity changes. The modelling was carried out for a chalk sample with an average porosity of 35%, with no shale content at a burial depth of 3100 metres. The synthetic time-lapse response in Figures 4.11 and 4.12 are generated using the calibrated rock physics model described in Chapter 2. Figure 4.11 shows that compaction resulting from water weakening and water flooding has the highest impact on P-wave and S-wave velocity changes. Smaller changes are observed due to dry compaction and pressure build-up. Unlike other clastic reservoirs, where pressure depletion in the oil leg is often not detectable above noise, in the Ekofisk field, the pressure depletion is coupled with porosity reduction; hence pressure depletion is often associated with a strong hardening signal. This, of course, depends on the pressure draw-down and the initial porosity of the rock. Figure 4.11 show how changes in P-wave velocity correlate well with compaction and saturation changes, whilst changes in S-wave velocity largely depend on pore pressure perturbation.

Similar observations are made for time-shifts, and the percentage change in P-impedance, compaction due to water weakening and sole water saturation changes plays a major role in increasing both of these attributes, as shown in Figure 4.12. Dry compaction, pressure build up and gas coming out of solutions present weaker signals in both the relative change of P-impedance and in time-shifts. In these plots, changes in time-shift and impedance due to gas saturation increases are not predicted to be as large as the changes induced by compaction or water saturation alone. This is because the initial gas saturation is already quite high (10%) in the LoFS surveys, therefore the non-linear portion of the property versus saturation behaviour curve is not accessed. Based on the prediction rom the fluid flow simulator, the gas saturation increase is minimal in the LoFS period, around a fraction of 0.3. These plots are useful to help in understanding the time-lapse signals measured from the LoFS surveys.

##### Figure 4.11: Showing the sensitivity of the percentage change of P-wave velocity and S-wave velocity to various production mechanisms that were modelled in isolation.

##### Figure 4.12: Shows the sensitivity of the percentage change of P-impedance and time-shifts to various production mechanisms that were modelled in isolation.

## 1.5 Derivation of a Proxy model via Synthetic Modelling

I have previously highlighted the rationale of a proxy model in various applications from reservoir characterization to seismic history matching. Most of this proxy model focuses on the relationship between the changes in dynamic properties and the difference in 4D seismic amplitude. Whilst the relation observed between dynamic properties and amplitude differences are encouraging, this particular approach is not the only one available, nor the most obvious choice for a thick, multi-cycle reservoir like Ekofisk. Another possibility is to relate dynamic changes to the relative change in impedances

∆IPIP×100

or other elastic properties which are interval properties. Impedances are usually modelled by a petro-elastic model to convert dynamic changes from the simulation model.

The relative change in impedance was chosen instead of seismic amplitudes, since it is also intuitive, and shows the relative strength of saturation or pressure changes in the 4D signals. The relative change in impedance is easier to interpret compared to the difference of impedance for baseline-monitor. The percentage change of impedance is confined within the range of -100% to+100%. It is more meaningful to look at this attribute than the difference in impedance, which is in the range of very large numbers. It is much more meaningful to quote that a 300psi change in pressure resulted in 1% change in impedance than 75000 m/s.kg/m^{3} change in impedance between baseline and monitor. Seismic amplitude is an interface or contrast property, whilst impedance is a layer or rock interval property; thus the latter is more easily related to dynamic properties. In the following sections, I will demonstrate the formulation of this proxy model using a synthetic model.

### 1.5.1 Synthetic model description

To analyse the individual impact of pressure, saturation and porosity reduction, the relative changes in impedances are modelled for a sector model of the Ekofisk field simulation model. This is carried out using the simulation-to-impedance procedure of Amini, MacBeth and Shams (2011), which from here on will be simply known as ‘sim2imp’. The sector model in Figure 4.13(a) has the same heterogeneity as the real data, and was simulated in a way that follows the production history of the actual dataset by limiting the number of wells to six injectors and six producers. Figure 4.13(b) illustrates the production evolution of field pressure, field oil, gas production rate, water cut and water injection rate. During the first 18 years of production, due to the nature of the oil and the initial pressure being close to bubble point, the gas came out of solution rapidly with poor pressure maintenance. In parallel to the real data, a water injection programme was initiated later to maintain pressure, resulting in an increase in oil production rate and water cut. Having simulated the dynamic changes, such as pressure, gas, water and oil saturation at different time steps, I then put these results into the sim2imp workflow to generate impedances.

In order to study the impact of compaction (dry compaction and water weakening) on pressure and saturation and, in turn, on the impedances, I ran the fluid flow simulator in two separate scenarios: (1) no compaction, and (2) compaction; this workflow is shown in Figure 4.13(c). In this approach, petroelastic parameters were firstly calibrated from the wireline logs (P- and S -wave velocities, density, water saturation and gamma logs) and fluid properties obtained directly from PVT data. The dry frame properties for the chalk were derived from laboratory rock-mechanics tests and stress-sensitivity curves coefficients were also taken from the laboratory and put into the stress-sensitivity model of MacBeth (2004). Fluid acoustic properties are calculated using black-oil PVT data in combination with the Batzle and Wang (1992) equations and are then mixed using a harmonic average. Calibrated chalk parameters are then used to perform fluid substitution. The petroelastic model used in the sim2imp procedure has been thoroughly described in Chapter 2.

Here, I demonstrate the results from my synthetic model for the individual contribution of pressure, saturation and compaction for case number 2, by including dry compaction and water weakening. Figure 4.14 (left panel) shows the dynamic changes such as gas, pressure and water by taking into account of compaction. Here, the term compaction encompasses physical changes such as dry compaction and water weakening. Figure 4.14 (right) shows the changes in P-impedance as percentages for the individual dynamic changes. Figure 4.14 (b – right, and c – right) depicts the effects of compaction due to pressure depletion and water invasion. The hardening signals in the relative change of impedance was amplified due to porosity reduction in those regions.

### 1.5.2 Reverse Engineering

Reverse engineering or backward engineering is, by definition, the process of extracting knowledge or design information and re-producing it based on that acquired information. In this case, I am trying to decompose the composite impedance from sim2imp to analyse its individual components and workings in detail, and recreate it using a simpler approximation. Modelling provides a way to examine the impact of pressure, gas, water saturation changes and geomechanical responses on the petroelastic parameters by independently isolating each of these controlling changes during the petroelastic modelling step, as shown in Figure 4.14. I can examine the effect of these different dynamic and geomechanical changes on the change of elastic properties, in the cases of no compaction, versus compaction. For simplicity the percentage change of any elastic properties are annotated in this chapter using the symbol ‘

δ

’:

δIP=∆IPIP×(100)

(4.1)

In the case of no compaction, I anticipated the composite impedance response might be decomposed as a sum of individual responses:

δIP∆P, ∆Sw, ∆Sg=δIP∆P,0,0+ δIP0,∆Sw,0+ δIP0,0,∆Sg

(4.2)

where

∆P, ∆Sw, ∆Sg

are changes in pressure and water and gas saturation respectively. By adding the independent

δIP

due to gas, water and pressure changes, the end product predicts the very same changes as if one ran running the full petro-elastic model. This is found to be accurate in the model across a wide range of geological and fluid conditions to within an error of 2% for P-impedance. This suggests that, in the case of no compaction, the system is additive for simultaneous pressure and saturation changes, and these changes honour the principle of superposition. This linearly additive behaviour forms the initial framework of the proposed proxy model equation.

The principle of superposition states that in all linear systems, the net response caused by two or more stimuli is the sum of the response caused by each stimulus individually. The relationship between the composite

δIP

from sim2imp is linearly related to the sum of the respective

δIP

computed from different dynamic changes. This is clearly demonstrated in Figure 4.15, with the cross-plotting of

δIP∆P,∆Sw,∆Sg

versus

δIP∆P,0,0+ δIP0,∆Sw,0+δIP0,0,∆Sg

. The goodness of fit is reflected in a high coefficient of determination,

R2

of 0.91.

R2

indicates how well a model obtained by linear regression fits the data. Therefore, the equation can be written in the form:

δIP=α∆P+b∆Sw+c∆Sg

(4.3)

The coefficients of

α, b

and

c

are obtained via linear regression of the individual dynamic changes and their respective impedance changes, while they are modelled in isolation of gas, water and pressure.

##### Figure 4.15: Showing the linear relationship between

##### IP∆P,∆Sw,∆Sg

and

##### δIP∆P,0,0+ δIP0,∆Sw,0+δIP0,0,∆Sg

for the case of no compaction.

The actual system contains compaction due to both dry compaction and water weakening, therefore, the same procedure is carried out by running the full simulator, but this time taking account of these geomechanical effects. In the sim2imp procedure, I then isolate the effect of each controlling factor, such as pressure, water saturation, gas saturation and porosity changes on the elastic parameters. In the water weakening system, the composite relative change in impedance can be similarly decomposed as:

δIP∆P’, ∆Sw’, ∆Sg’,∆φ=δIP∆P’,0,0,∆φ+ δIP0,∆Sw’,0,∆φ+ δIP0,0,∆Sg’,∆φ

(4.4)

where

∆φ

represents the change in porosity; this extra term is attributed to compaction.

Additionally, the impact of compaction on the changes of gas reflected as relative changes in P-impedance inside the reservoir can be studied by taking the difference of the functions derived from the compaction case and those from the no compaction case for gas:

δIP0,0,∆Sg’,∆φ-δIP0,0,∆Sg,0

. The difference between

δIP0,0,∆Sg’,∆φ-δIP0,0,∆Sg,0

in Figure 4.16(a) shows compaction results in less gas liberation and in turn an increase in the relative change of P-impedance. Likewise, for the impact of pressure on the relative change in P-impedance as a function of compaction, in Figure 4.16(b):

δIP∆P’,0,0,∆φ -δIP∆P,0,0,0.

Thus, as a result of pore collapse and higher compressibility, there is an increase in pore pressure and an increase in relative P-impedance. Moreover,

δIP(0,∆Sw’,0,∆φ )-δIP0,∆Sw,0,0

in Figure 4.16(c) shows that the intrusion of water results in further pore collapse and reduction in permeability, which causes a greater overall relative change of P-impedance. Figure 4.16(d) shows the changes in porosity due to compaction. The change in porosity

,∆φ

is given as the difference between initial porosity and the updated porosity (in the monitor survey). In comparing the influence of compaction on pressure, water and gas, I found the influence on pressure and water is much greater, with approximately 3% increase in relative impedance, while for gas it was slightly below 0.2%.

##### Figure 4.16: Cross-sections for the sector model, showing the

##### δIP

between the cases of compaction and no compaction for (a) gas, (b) pressure (c) water saturation and (d) porosity change.

To assess if the composite impedance,

δIP∆P’, ∆Sw’, ∆Sg’,∆φ

is indeed linearly correlated to the sum of the individual impedances from dynamic and porosity changes,

δIP∆P’,0,0,∆φ+ δIP0,∆Sw’,0,∆φ+ δIP0,0,∆Sg’,∆φ

, both terms are cross-plotted in Figure 4.17. The best fit, linear regression yields a goodness of fit,

R2

of 0.85, which shows that these terms are linearly related.

##### Figure 4.17: Showing the best fit, linear relationship between

##### δIP∆P’, ∆Sw’, ∆Sg’,∆φ

and

##### δIP∆P’,0,0,∆φ+ δIP0,∆Sw’,0,∆φ+ δIP0,0,∆Sg’,∆φ

In addition, I assessed the linearity behaviour between compaction versus the change in impedance resulting from porosity reduction. This linear relation is depicted in Figure 4.18, by cross-plotting the effect of compaction on

δIP(∆P’,0,0,∆φ)

versus the changes in porosity,

∆φ

. The linear proxy function for the compaction case can be written as:

δIP=α∆P’+b∆φ+c∆Sw’+d∆Sg’

(4.5)

##### Figure 4.18: Cross-plotting of the change in porosity

##### ∆φ

versus the relative change in impedance due to compaction

##### δIP(∆P’,0,0,∆φ)

also yields a linear relationship.

In the case of compaction, by adding the independent relative impedance changes due to gas, water and pressure changes, as shown in Figure 4.19(b), the end product predicts the same change as through the full petro-elastic model, as illustrated in Figure 4.19(a). This shows that the principle of superposition is also applicable in the case of compaction, and that the approximation of the petro-elastic model as a single equation is valid for the modelling of relative change in P-impedance.

##### Figure 4.19: Comparison of the cross-sections for the composite impedance difference in percentages between baseline and monitor from sim2imp to the sum of individual difference in impedance, in percentages.

### 1.5.3 Linearisation of compaction curves

Apart from providing a faster way to compute the relative changes in elastic properties, the other objective of this proxy model is to estimate porosity reduction efficiently and intuitively. This is achieved by incorporating rock mechanics laboratory data as constraints into the proxy model equation. The crucial parameters that describe porosity reduction can be computed via a set of compaction curves for dry compaction and water weakening. These curves are similar to those used in the fluid flow simulation model. The changes in porosity can be described as the difference between porosity during the initial and subsequent time steps:

∆φ=φi-φupdated

(4.6)

Here,

φi

represents the initial porosity and

φupdated

represents the porosity after simulation conditions, whereby

φupdated<φi

, as the compaction process is irreversible.

φupdated

can also be expanded as the partial differentiation of porosity to changes in pressure:

∆φ=φi-∂φ∂∆P∙∆P+φi

(4.7)

∆φ=-∂φ∂∆P∙∆P

(4.8)

The compaction curves for both dry compaction and water weakening mechanisms provided by the field operator are displayed in Figure 4.20(a-left) and 4.20(a-right). These curves were calibrated using well data to interpolate the compaction gradient for each of the initial porosity classes. In work from Sylte et al. (1999) and Janssen, Smith and Byerley (2006) radioactive marker bullets were instrumented in wells and repeatedly logged to monitor strain and displacement in relation to reservoir compaction, and these data were used to calibrate the compaction model. I first linearized the compaction curve for each porosity curve in both dry compaction and water weakening behaviours, essentially replacing the gradient for each initial porosity member with a constant. Using a stepwise algorithm, I replaced

∂φ∂∆P

by two numerical functions called

Fp

and

Fww

, each describing the gradient for initial porosity based on the amount of change in pressure; this captures the reduction of porosity for both dry compaction and water weakening. These functions are depicted in Figure 4.20(b).

Fp=-∂φ|∆P∂∆P

(4.9)

Fww=-∂φ|∆sw∂∆P

(4.10)

Having replaced the porosity reduction terms in Equation (4.5) with the compaction functions,

Fp,

and

FWW

, the modified equation is given as:

δ∆IP=α∆P’+b(Fp+ FWW)∆P’+c∆Sw’+d∆Sg’

(4.11)

The

Fww

case occurs when the maximum attainable water saturation (

Swmax

) is reached. The model was first put forward by Sylte et al. (1999), showing no additional compaction for

Sw

> 0.325 when water invasion occurs. Therefore, even when the water saturation exceeds the maximum attainable water saturation (

Sw

=0.33), there will be no additional compaction (Das et al. 2016). However, in many cases, the increase in water saturation does not exceed this value, and hence a new function is introduced to avoid overestimation of the relative change in impedance due to maximum water weakening. This pseudo function

F̂ww

manipulates the correct amount of compaction of the two end members

Fww

and

Fp

by weighting each one with the ratio of

∆Sw+SwiSwmax

; this function is shown in Figure 4.20(b). The pseudo-function is described as:

F̂ww=∆Sw+SwiSwmaxFww-Fp

(4.12)

The equation can be re-written by including the pseudo function

F̂ww

to replace the end member case of

Fww

. The proxy model equation is now given as:

δIP=α∆P’+bFp+∆Sw+SwiSwmaxFww-Fp∆P’+c∆Sw’+d∆Sg’

(4.13)

Since the initial water saturation,

Swi

in the Ekofisk field is low, at 0.05, it can be treated as negligible and the equation is further simplified as:

δIP=α∆P’+bFp+∆SwSwmaxFww-Fp∆P’+c∆Sw’+d∆Sg’

(4.14)

In summary, the proposed proxy model is derived analytically through modelling and reverse engineering and a comprehensive synthetic modelling exercise. To predict the relative changes in impedance, the proxy model requires a set of calibrated coefficients (

a

,

b

, c and

d

) from the petro-elastic model, changes in pressure (

∆P

) and saturation (

∆Sw

), and compaction curves from the laboratory, expressed as functions (

Fww

and

Fp

) of initial porosity

φi

. The procedure used to construct the proxy model for the relative change in P-impedance,

δIP

is also carried out for S-impedance,

δIS

, P-wave velocity,

δVP

, S-wave velocity,

δVS

, and density,

δρ

. A different proxy model equation for each of the elastic parameters can be derived separately, where the general form is similar but coefficients

a, b, c

and

d

have different magnitudes.

### 1.5.4 Validating the proxy model

One way to cross-check the validity of the proxy model, is to compare the porosity reduction estimated from the fluid flow simulator with that calculated from the proxy model. The input parameters required in the proxy model are the changes in pressure and water saturation, and the initial porosity; these are readily available as outputs from the fluid-flow simulator. Figure 4.21 shows a comparison between the results from the simulator (left) and the linearised compaction functions (right). In general, I obtained a good agreement between both estimations, with marginal discrepancies. The unevenness of the porosity change generated from the fluid simulator is caused by rock typing in the model. As previously mentioned in Section 4.2.1, the different rock types are assigned with different relative permeability and rock compaction information although they can have the same porosity and permeability ranges. The compaction functions on the other hand, are irrespective of rock typing and the compaction gradients employed are only a function of the initial porosity.

##### Figure 4.21: (Left) porosity reduction simulated from ECLIPSE compared to (right) porosity reduction calculated from

##### FP

and

##### FWW

Another way to gauge the accuracy of the estimation from the proxy model is to calculate the percentage error between the percentage change in P-impedance from the sim2imp procedure and the proxy model for the individual components of gas, water and pressure. The percentage error is described as:

δIPsim2imp-δIPproxyδIPsim2imp×(100)

(4.15)

The percentage error between the full sim2imp procedure and the proxy model is generally small, with errors less than ± 5%, with occasional high points at ± 1.5%. The percentage error is shown in Figure 4.22 (a, b, and c) for the respective changes in P-impedance as a function of gas and water saturation and of pressure change.

##### Figure 4.22: Showing percentage error of the percentage change of P-impedance estimated for (a) gas saturation, (b) water saturation and (c) pressure change between sim2imp and proxy model

## 1.6 Solving the Inversion Problem with a Proxy Model

I will next demonstrate how the proxy model is applied in inversion. Modelling and inversion are closely related: in modelling we seek to reproduce an observation (or measurement) by perturbing parameters that are somehow related to such observations. Inversion uses a series of observed measurements to calculate those parameters we are interested in, by calculating a series of predictions from an established model and comparing them to the observations (Menke, 1989). In both cases, the core of the process is to establish a model which relates the measurements with the parameters we wish to estimate and vice versa. In this section, I would like to demonstrate the following propositions:

- The inversion scheme using the proposed proxy model produces realistic results for the synthetic data.
- Constraints and prior information are crucial in narrowing the search space and subsequently provide us with better estimates of the model parameters.

### 1.6.1 The inverse problem and the optimisation solution

The inverse problem, as described by Menke (1989), can be mathematically defined as an integral function that relates the measurements

di

to the parameters we wish to estimate, (x), through ( ), which relates the two:

di=∫G(x)m(x)dx

(4.16)

where

d

is the data observation vector and

m

is given as the parameter vector and

G

is the data kernel matrix. Here, I attempt to solve the inverse problem on a sample by sample or cell by cell basis in the simulation model, the inverse equation can be written as:

di=∑j=1Gijmj

(4.17)

The vector

di

represents the synthetic elastic properties: the percentage changes in P-impedance, S-impedance, P-wave velocity, S-wave velocity and density generated using the sim2imp procedure. The parameter vector,

mj

, contains the pressure and saturation changes to be estimated and the data kernel matrix,

Gij

, is the forward operator relating the two. Given the proposed proxy model equation takes a non-linear form, this inherently makes the inverse problem particularly difficult. There are two reasons: firstly that nonlinear error propagation is a difficult problem, and secondly, it introduces non-uniqueness in the solutions. If the forward problem is linear, and if a *L*2-norm is used for the calculation of the misfit, the misfit function has a parabolic dependence on the model parameters, and therefore the misfit function has a single minimum. Any type of descent method will lead to this unique minimum. When the forward problem is non-linear, the misfit function can have multiple minima. The problem with these local minima is that search methods for the global minimum may misidentify a local minimum as the global minimum. In that case the estimated model is not the model that gives the best data fit.

The idea of an optimisation is to achieve the best possible result in acceptable conditions. Since the problem is described as non-linear, I cannot solve this via the linear least squares solution and need to tackle this deterministically, using other numerical optimization solutions. Typically, one seeks to recover the model parameter,

m,

based on observations,

d

, where both are related by a forward modelling operator

G

, as discussed earlier. The problem is to find *m* such that the misfit is less than a certain tolerance:

minimizeG(m)-d

(4.18)

I start off by assuming that a local optimal solution exists for this problem, and a local search algorithm will suffice for my cause. I will talk more about using global algorithms in the subsequent section; this sort of algorithm is computationally expensive and is capable of looking at objective function values more exhaustively in the search space. The traditional local algorithms are subdivided into two main classes: direct (search) and gradient-based methods (Reklaitis et al., 1983, Deb,1995, 1998). In direct methods, only the objective function and constraints are used to guide the search strategy. The gradient-based methods (e.g. linearized inversion) use the first and/or second-order derivatives of the objective function and/or constraints to guide the search process. These algorithms can converge quickly to the solution. For quasi-linear problems, these algorithms are a good selection. I will be employing the Trust-Region Reflective algorithm available from the MATLAB optimisation toolbox, which is a subclass of the gradient-based methods. The Trust-Region Reflective algorithm used in the optimization is described in Coleman and Li (1994, 1996). Each iteration involves the approximate solution of a large linear system using the method of preconditioned conjugate gradients (PCG). A typical simple definition of the objective function is through the sum of squared differences:

f=∑k=1Nd(dkobs-dkmod)2

(4.19)

where the sum is taken over all available observed data,

dkobs

, and

Nd

is the total number of these data. The objective function of this problem is solved cell-by-cell, it can be written in the same format as:

minimize f∆P,∆Sw,∆Sg

where

f∆P,∆Sw,∆Sg=δIP-GIp∆P,∆Sw,∆Sg2+δIS-GIs∆P,∆Sw,∆Sg2+δVP-GVp∆P,∆Sw,∆Sg2+δVS-GVs∆P,∆Sw,∆Sg2+δρ-Gρ∆P,∆Sw,∆Sg2

(4.20)

G

is the forward operator and represents the proxy model equation, which is unique to generate different elastic properties such that the coefficients in the proxy model for P-impedance will be different from S-impedance, P-wave velocity, S-wave velocity and density. The inversion of various data sets was done jointly. Joint inversion is carried out in this case to produce mutually consistent estimates of the unknown parameters. As demonstrated in Equation 4.20, one objective function is to be optimized from the summation of individual objective functions representing various data sets. In this synthetic example, given that the data is generated and not measured from the field, the quality of the data is consistent. However, for the actual field data, the quality of the different types of measurement may differ, such as for time-shifts, amplitude and the inverted elastic properties. The component objective functions should therefore be multiplied by weight factors giving them the correct contribution for determining the model parameters (Drahos, 2008). This is discussed further in the next chapter when using actual field data. As previous researchers (Floricich, 2006; Alvarez, 2014) observed, in obtaining reliable estimations of

∆P

and

∆Sw

, especially with the highly correlated nature of the 4D signals, incorporating additional constraints into the inversion scheme is helpful. There will be several possible combinations of Δ and Δ which will produce the same change in petro-elastic properties, for example, both a decrease in pore pressure and an increase in water saturation will generate a hardening signal. The implication of this is that the inversion problem is ill-posed and hence it is required to provide constraints in order to find a solution that falls within the expected ranges of change.

### 1.6.2 Constraints

To find a unique solution, one must add some constraints in the solution space. For example, in order to aid the convergence of the optimisation, a good initial guess or starting point in solving the optimisation problem is helpful. In this case, I used the mean value from the simulation model prediction for gas, water and pressure changes (

∆P̅,∆SW, ̅ ∆Sg̅

) as the initial starting point.

Another way to ensure a better feasible solution is by setting boundaries for the parameter estimates. If one knows the bounds on the location of an optimum, one can obtain faster and more reliable solutions by explicitly including these bounds in the problem formulation. For example, the change in water and gas saturation cannot exceed 0 and 1, respectively, with further constraints in certain compartments of the reservoir, such that no water saturation changes should be expected in the water leg, since the water saturation is already unity or at its maximum. Moreover, it is unlikely for other types of reservoir fluid to replace water in the aquifer due to density variation, where water is denser than oil and gas. On rare occasions, the injectors could push oil into the aquifers, however this is not the norm and will not be considered in this modelling. In most parts of the reservoir, the water saturation change is constrained to:

0≤∆Sw≤(1-Sor-Swirr≈0.675)

. In the Ekofisk field, the irreducible water saturation,

Swirr

is given as 0.05 and the residual oil saturation,

Sor

has a value of 0.275. Below the oil-water contact, the change in water saturation is

∆Sw=0

.

Based on the prediction of the simulation model for gas changes, the maximum value peaks at 0.6, and no gas was expected to go back into solution, due to a significant overall pressure depletion in the reservoir. Therefore, the expected gas saturation change is limited to

0≤∆Sg≤0.6

. On the other hand, the pore pressure change in the reservoir can also be constrained by studying the production history of the field. Based on the minimum and maximum estimates of the change in pressure predicted from the simulation model, the bounds set on pore pressure change is between:

-20≤∆P (bar)≤0

. These constraints can be incorporated into the optimisation problem in the form of an equation such as

Hm≥h

:

100-1000100-1000100-1∆Sw∆Sg∆P≥00.67500.6-200

(4.21)

Since I know the location of the original oil-water contact prior to production, the constraints in the water leg are given as the equations below. Below the oil-water contact, the water saturation will always be 1, and no change is assumed to occur; the model parameter for

∆Sw

equals to 0. Below the oil water contact, the optimisation problem is solved subject to these conditions:

-20≤∆P (bar)≤0

(4.22)

∆Sw=0

,

∆Sg=0

(4.23)

## 1.7 Results and Discussion

The non-linear inversion solution is a non-unique one, given there are a multitude of models that explain the data equally well. My approach here is to reconstruct an estimated model (proxy model) that still captures the full physics as the true model, for the use of forward modelling and inversion. The proxy model is capable of forward modelling a series of different elastic properties with good agreement to the data estimated from the true model, which in this case is represented by the sim2imp procedure. Of course, one can also argue that a different rock physics model will give a different set of data. In this case, the true model is assumed as the rock physics model calibrated with well measurements in the sim2imp procedure, which is used as a benchmark for estimation of elastic properties. In this section, I will carry out two separate inversions, highlighted in two separate routes as ‘1’ and ‘2’ in Figure 4.23. The first route involves using the data from the proxy model to invert for the model parameters. The second route uses data estimated from the full physics model but inverted back to the respective model parameters using the proxy model.

Figure 4.24(a) shows the comparison between the dynamic changes predicted from the simulation model and Figure 4.24(b) shows the inverted dynamic changes using the proxy model from data generated from the same proxy model in route 1. Route 1 is also known as ‘inverse crime’ (Wirgin, 2004), when the same (nearly same) theoretical ingredients are employed to synthesize as well as to invert data in an inverse problem. In general, route 1 managed to reproduce the pressure and saturation changes, even in places where these changes overlap (close to the crest of the reservoir, where gas accumulates and pore pressure reduces). The model data are uniquely resolved with very small residual errors. This is expected, as the same data kernel (G’) is used in both forward and backward modelling. This procedure is useful to investigate the stability of the inversion solution in the presence of noise and to assess the effectiveness of the proxy model in an underdetermined problem. To truly assess the robustness of the proxy model, the inversion is carried out using data generated from the full physics model, as described in route 2.

##### Figure 4.23: The workflow comprising the forward modelling and inversion routes using the full physics model (sim2imp) versus the proxy model.

##### Figure 4.24: Cross-sections showing (left) dynamic changes generated from simulation model compared to (right) inversion results using proxy model for (a) gas, (b) water saturation and (c) pressure change. The input data for the inversion are generated from the proxy model equation.

Inversion results from the second route are displayed in Figure 4.25(a, b, and c). In general, the inversion results shows greater discrepancy between the inverted model parameters and with those simulated from ECLIPSE, as compared to the results shown in Figure 4.24(a, b, and c). This is particularly apparent in places where there are overlaps between pressure and saturation changes, such as the gas cap, and around the injector well I4. However, on a positive note, overall the polarities of the inverted dynamic changes are comparable to those simulated from ECLIPSE, which is essential for interpretation. These discrepancies are expected, since a different data kernel is used in forward and backward modelling.

The residuals and the convergence criteria are assessed in this section. The residual error is the value of the objective function at solution, which can be calculated for the independent observed data (

δIP, δIS, δVP,δVS, δρ

). On the other hand, the exit criteria show the reason for the solver to stop the calculation. These exit criteria are unique in the MATLAB optimisation package. When the optimisation solver completes, it sets an exit criterion. The exit criterion is an integer that describes the reason the solver has halted its iterations. In general, positive criteria corresponds to successful solutions and negative ones are not. The zero exit criterion represents a solution that is computed by exceeding an iteration limit or limit on the number of function evaluations (MATLAB documentation). The exit criteria are colour-coded from -2 to 4 and explained in Table 4.3.

#### Table 4.3: The exit criteria and their description for the non-linear optimisation solver (MATLAB documentation).

Figure 4.26(a) to (e) shows the residual error (in route 2) for each of the predicted elastic properties at the objective function when a solution is found, and Figure 4.26(f) represents the exit criteria for each cell in the model. I observed that the residual errors for the different elastic properties are generally small, with higher residuals in composite elastic properties such as P and S-impedance. Because these result from the multiplication between P, S-wave velocity and density, the residual error is also amplified in these elastic properties. The residual errors are expected because the coefficients derived for the proxy model are calculated by linear regression between elastic properties from the sim2imp and simulation estimated dynamic properties, which might not be able to capture all the scatter in this approach. However, this regression is considered a robust calibration, as in the derivation of the proxy model and its coefficients, all data points were used in the regression exercise.

The exit criterion show mostly ‘1’, which informs us that the solver stopped because it has converged to a solution. In some cells, where the exit criterion is numbered as ‘3’, this shows the change in the residual norm is very small; hence the solution is accepted. Given this is a non-linear problem, this could suggest that there could be multiple local minima. Thus, a different solution can be applied, such that one could make inferences from a range of models and search in many directions. This brings us to the next section, where I suggest an alternative solution, although the current solution shown in Figure 4.26 has already yielded minor discrepancies compared to the true model parameters.

### 1.7.1 Generation of populations that fit the data

An alternative approach is to compute the misfit for a very large class of models and to use the data fit, possibly in combination with Bayesian statistics, to make inferences about the range of models that explain the data in a certain sense (e.g. Mosegaard and Tarantola, 1995, Gouveia and Scales, 1997, 1998, Mosegaard 1998). Obviously, this approach requires a numerical approach to create such ensembles, but if the forward modelling can be easily computed, like the proxy model, computation time can be greatly reduced. An important concept in the generation of ensembles of models is the randomness in the search method that one employs. This stochastic approach will be applied to the actual data in Chapter 5.

## 1.8 Summary

This chapter first addressed the need for a proxy model to invert for pressure and saturation changes in the impedance domain. Pressure and saturation changes are more easily decomposed from the elastic properties rather than from seismic amplitude differences. The aim of the remainder of the chapter was to show the derivation and the validation of the proxy model on a set of synthetic data.

The proxy model derived in this chapter is set to achieve the aforementioned goals, which are to provide a fast-track method to both forward model elastic properties and invert for pressure and saturation changes in a compacting chalk reservoir. The analysis reported in this chapter has validated the competence of the proxy model approach for both forward and backward modelling, by comparing the discrepancy in estimated elastic properties with the true model, which is the sim2imp procedure. Marginal differences were observed between the predicted elastic properties between the proxy model and those calculated by sim2imp. Similarly, the porosity reduction calculated using the proxy model was also comparable to those predicted from the ECLIPSE simulation model. This confirms that the compaction behaviours can be simplified into analytical equations and be embedded in the proxy model.

The inversion problem is cast as an optimisation problem, where for each cell, multiple forward modelling is allowed, to minimize the misfit between the predicted data and the observed data. Crucially, by including constraints such as bounding values and statistical information from the simulation model, the optimisation process to solve for pressure and saturation changes is expedited. I have kept the use of prior information to the minimum. In practice, the selection of bounding constraints and prior information depends on data availability. For the actual data I will also explore possible a priori information from a well history-matched simulation model and engineering concepts. The inversion using the proxy model equation shows promise, and the inverted solutions are also considered stable. The residual errors are less than ±2% for all estimated data.