Scale and Process Dependent Model Errors in Seismic History Matching

Résumé — Erreurs résultant des effets d'échelle et de la simplification des processus dans le calage de modèles à partir de données sismiques — Les données sismiques fournies par la méthode timelapse (4D) renvoient des informations spatiales et dynamiques sur les modifications des propriétés du fluide du réservoir et peuvent être utilisées pour imposer des modèles de simulation de flux, améliorant ainsi la confiance en la caractérisation du réservoir et son comportement prévu.


INTRODUCTION
Reservoir management requires tools such as simulation models to predict asset behaviour.History matching is often employed to alter these models so that they compare favourably to observed well rates and pressures.This well information is obtained at discrete locations, and thus lacks the areal coverage necessary to accurately constrain dynamic reservoir parameters such as permeability and the location and effect of faults.Time-lapse seismic captures pressure and saturation changes by their effect on seismic impedances or equivalent attributes.This information, missing in conventional approaches, can be included in seismic history matching improving estimates of the reservoir model parameters.
We have developed such a method (Fig. 1) that is automated, based on a multiple model, quasi-global stochastic approach and has been applied to the Schiehallion field, located to the west of the Shetland Islands on the United Kingdom Continental Shelf (Stephen et al., 2005b;Stephen and MacBeth, 2006a;Stephen and MacBeth, 2006b;Stephen and MacBeth 2006c;Stephen, 2006;Stephen et al., 2006).Our method improves on the classical workflow, where the engineer manually adjusts parameters in the simulation model.It also improves on gradient-based methods, such as Steepest Descent, Gauss-Newton and Levenberg-Marquardt algorithms (e.g.Lépine et al., 1999;Dong and Oliver;2005;Gosselin et al., 2003;Mezghani et al., 2004), which are good at finding local misfit minima but can fail to find the global minimum.Our approach is also faster than stochastic methods such as genetic algorithms and simulated annealing, which often require more simulations and may have slower convergence rates.Similarly our method is more efficient at searching the parameter space than other probabilistic methods such as the Gradual Deformation (Hu, 2000) or Regional Probability Perturbation (Hoffman and Caers, 2005) methods.Finally, multiple models are generated enabling posterior uncertainty analysis in a Bayesian framework (e.g.Stephen and MacBeth, 2006a-c).
Our approach is based upon the prediction of well behaviour and also changes in seismic attributes during the production process.For the purpose of this paper, we shall refer to this prediction process as 'the model' which itself consists of several modelling steps.The first of these involves the generation of models that represent the geological architecture of the reservoir.These are commonly referred to as geo-models in the literature.Following any necessary upscaling, geomodels may be used as input for flow simulations.This second modelling step can be used to predict fluid changes under production conditions and both geo-modelling and flow simulation are standard practice in the oil and gas industry.A third modelling step uses the results of the simulations as input into a petro-elastic transform to predict changes in the petro-elastic properties of the reservoir.We compare these to equivalent properties from the observed seismic  seismic data offers spatial and dynamic information about changes in reservoir fluid properties and can be used to constrain flow simulation models thereby improving confidence in the reservoir characterisation and its predicted behaviour.To address this, we have developed a method of quantitatively integrating 4D seismic data in an automated history matching workflow.Appropriately parameterised flow simulations are converted to predictions of 4D signatures by a petro-elastic transform and suitable rescaling before a misfit is calculated by comparison to observed data.Model parameters are then updated using a quasi-global stochastic inversion method.This process is affected by scale and process dependent model errors.Flow simulations are often created such that computer resources are optimised and some level of accuracy is sacrificed.To speed up simulations, some form of upscaling is required to capture two-phase flow properties such as relative permeability but also to represent geological heterogeneity.The upscaling may be over-simplified or ignored.In addition, simplifications to the flow processes may be made, for example by using streamline methods.Finally, the petro-elastic transform contributes to the model errors due to assumptions about saturation distributions and cross-scaling is required because modelled and observed seismic are obtained for different volumes.We present an analysis of the above model errors that occur using a synthetic geo-model based on a North Sea reservoir.We show that the model error depends on the rock physics parameters as well as the underlying geo-model.When the 4D signature is dominated by pressure effects, the model error is negligible in our case.We describe how the model error affects the history matching process due to biasing.The latter results in a best set of model parameters which may be different from that obtained by upscaling while the uncertainty estimator is also changed.We compare the effect of the model error to other errors such as observed data errors.Finally, we describe how the model error is addressed in the misfit calculation to improve the history matching process and reduce the biasing effect.

Abstract -Scale and Process Dependent Model Errors in Seismic History
although they could be used as input into seismic modelling (e.g.MacBeth, et al. 2005).Our model therefore combines geo-modelling, flow simulation and the petro-elastic transform with their associated errors.
To date, we have ignored the influence of model errors in the history matching process.These can be important, however, as they can affect the convergence rate, lead us to incorrect models or alter the posterior uncertainty estimate.The errors arise because we make simplifications in the simulation process as discussed in Christie et al. (2005).One such simplification occurs in the physical processes that we model.For example the complex pore scale behaviour is represented by continuum equations and fluid interactions are captured by relative permeability and capillary pressure functions.We also convert those equations to a simplified form by discritisation via finite difference methods.Our equations become less accurate the more coarsely we discretise but we also change our representation of the finer detail in both parametric form and in the physics we hope to capture.This applies to both flow simulation and the conversion of simulation results into predicted seismic (e.g.Stephen et al., 2005).
The purpose of this paper is to calibrate the scale dependence of the model errors to determine how we may remove or reduce their effect.We do so by first assuming that the model is accurate at some fine scale and derive synthetic predictions of the reservoir behaviour.Our model is analogous to the Schiehallion field model that we studied previously.We then generate appropriate coarse scale models and attempt to modify these using our history matching method so that predictions compare to those made at the fine scale.We show that the model error affects our ability to history match to a degree that is equivalent to errors that may be present in the data.We can remove the impact of the model errors, however, with effective calibration in the misfit calculation similar to O'Sullivan (2004).

METHOD
Our approach may be summarised as two separate steps common to all iterative history matching methods.First we take a set of models where the p-th model is parameterised by the vector m p and calculate a misfit, J(m p ) for each model.Then we use a method of choosing a new set of m p based on the existing misfits, in this case the Neighbourhood Algorithm (Sambridge, 1999a).We then iterate through the two steps until we converge on a set of good models.
The first step is quite complex and involves several modelling steps and so we have broken up the description appropriately (Fig. 1).The Neighbourhood Algorithm is independent of these steps and we give a brief summary here, although further details can be found elsewhere (Sambridge, 1999a).We begin the process creating an initial ensemble of models by sampling the parameter space either randomly or deterministically.

Generating Multiple Simulations
For the p-th parameter vector, m p , a new geo-model is derived by modifying a base case, created using sequential Gaussian simulation using the GSLIB package (Deutsch and Journel, 1992).In this paper we only vary the permeabilities in the geo-model using the Pilot Point Method with Kriging (e.g. de Marseily et al., 1984).The vector, m p , is therefore a set of permeability multipliers applied to the pilot points, specific cells in the geo-model.Kriging is used as a form of interpolation to change permeability multipliers between pilot points.The variogram used in the Kriging has a range equivalent to the separation of the pilot points.In the Schiehallion field study (e.g.Stephen and MacBeth., 2006;Stephen et al., in Press) we modified additional parameters including the Net:Gross distribution, to match pre-production acoustic impedance, as well as fault transmissibilities and the petroelastic transform parameters when matching dynamic data.

Forward Modeling: from Flow Simulation Results to Seismic Impedances
The observed seismic is compared to the predicted impedance attributes, which are obtained using a petro-elastic transform (e.g.MacBeth et al., 2005) and fluid saturations and pressures from the simulations.For each simulation cell, we calculate the dry bulk modulus using (MacBeth, 2004): (1) where the superscript r identifies rock type (sand or shale), the parameters κ inf , E κ and P κ are determined from lab measurements or by history matching (Stephen and MacBeth, 2006c), and represents the dry bulk modulus at Standard Temperature and Pressure, the excess compliance present in the rock as a result of geological or mechanical processes, and the stress sensitivity respectively.P eff is the difference between the overburden pressure and the pore pressure.Here we assume that the effective stress equals the differential.The shear modulus, μ r , has the same form with equivalent parameters.
We then use Gassmann's equation (1951) to get the saturated bulk modulus for each lithofacies within a cell: (2) where φ is the porosity, κ gr is the bulk modulus of the mineral, α = (1 -κ r dry /κ gr ).κ f is the fluid bulk modulus given by the saturation weighted harmonic average of the individual phase bulk moduli (e.g.Domenico, 1976) from: (3) where S w , S o and S g are the water, oil and gas saturations respectively and κ w , κ o and κ g are the water, oil and gas moduli respectively and obtained using Batzle and Wang (1992).
The P-wave moduli for sand and shale are obtained from M r = κ r sat + 4 μ r /3 (shale is assumed to consist of dry frame only and the shear modulus is unaffected by saturation) and the value for each cell, M cell , is obtained from the harmonic mean of the sand and shale values, weighted by the respective fractional volumes via the Net:Gross.This is valid for vertical propagation in a layered medium (Backus, 1962).Using this, the impedance for a column of cells in the simulation model is calculated: where ρ is the bulk density of the cell obtained by averaging the densities of the rock frame and the fluid densities.The brackets, "< >", indicate a vertical volume weighted average over the reservoir interval.This approach is valid for reservoir beds that are less than one tenth of the seismic wavelength thick and units of around one quarter wavelength thick.By averaging over the reservoir interval in this way, we get a map of predicted impedance with areal resolution at the scale of the simulation grid.
For the Schiehallion field study, we calculated 4D attributes from the observed data by integrating over the reservoir interval.In the synthetic study presented in this paper, we synthesise the observed data by using impedances calculated from a fine grid flow simulation model with cells Comparison of the coarse and fine grids used in our study.Thick grey lines indicate the coarse cells while large circles show the location at which the impedances are predicted.Equation 5 is used to interpolate the impedances to obtain values at the small black circles, i.e.where the observed seismic would normally be measured.Broken and solid arrows indicate the principal directions of the coarse (simulation) and fine (seismic) grids respectively.

Simulation (J)
Cross lines (j) Simulation (I) measuring 20 × 20 × 1.2 m.During history matching we run simulations on a coarse scale grid measuring 100 × 100 × 6 m.The predicted impedance is obtained at a coarser scale compared to what we observe, therefore.To enable comparison we must either upscale the observed data or downscale the predicted impedances.We prefer to keep our observed data as intact as possible and interpolate the predicted seismic from (Fig. 2): (5) where indices I and J are the x and y indices on the simulation grid and i and j are equivalent for the finer seismic grid, is obtained from Equation 4and w ijIJ = exp(-β(|r IJr ij |) and r is the position vector of the centre of the cell.We found that β = 0.05 m -1 gave the best results in our field study, minimising the representivity error.
The same process was applied in the Schiehallion field study (Stephen et al., 2005 onwards) except that there, the observed and predicted seismic data grids were aligned at 45 degrees to each other whereas here they are parallel.In that study, the bin size of the observed seismic was 12.5 m.

Evaluate Misfit
A single misfit objective function is obtained for each model incorporating a comparison between observed and predicted production and seismic data.For each variable being compared we use the following equation (Tarantola, 1987): where the x i is the i-th data vector being compared, with superscript o or m for observed or modelled respectively, and is the covariance matrix of x i and is the sum of covariance matrices of the observed and model errors, and respectively.
The total misfit per model, J, is then the sum of misfits for each variable: (7)

Error Analysis
We can split the observed and modelled variables (e.g.Glimm et al, 2004) where x i t is the true value of the variable, and ε i is the error and superscripts o, m and p correspond to the error of the observed data, the model error and parameter error.We aim to reduce the parameter error to zero by history matching.
Seismic attributes contain errors from a number of sources including repeatability of the acquisition components, crossequalisation, tuning, wave interferences, and time integration which are difficult to assess.In the Schiehallion field (Stephen et al., 2005b, etc.), we determined the error by calculating the covariance matrix, directly having calculated ε i o using a Weiner band pass filter (Press et al., 1998) to identify the signal.In this study we use the noise information from that field study to perturb the synthetic impedance so that we can analyse the relative impact of the errors.We also combine data conversion errors such that they are not strictly related to measurements.
In general in reservoir engineering, we perform simulations using the finite difference approximation of the flow equations together with a discretised rock volume and this approach inevitably introduces some errors.These exist no matter what scale we simulate at although we generally expect that they are reduced as we refine the simulation grid.The simulation model is often based on a finer grid geomodel which has been upgridded and upscaled.Parametric errors may be generated if we upscale the original geo-model incorrectly so that we fail to represent sub-grid flow properly (e.g.Christie et al., 2001).We also get numerical dispersion (spreading of the simulated flood front) if we use relative permeability and capillary pressure curves measured at the wrong scale (Kyte and Berry, 1975) or, alternatively, other model errors appear if we upscale these functions incorrectly (Barker and Dupouy, 1999).
Another source of error may arise when we calculate the impedances.Domenico, (1976) observed the saturation relationship in Equation 3. Researchers have since pointed to the scale dependence of this equation (e.g.Knight et al., 1996).If we apply it to large scales (i.e.large grid cells) then the seismic induced pressure wave is not fully equilibrated over the sample volume.We can estimate the maximum cell size from the critical length, λ d , (e.g.Mavko and Mukerji, 1998) derived as: where k is the permeability, F is the seismic frequency, μ, is the viscosity of phase, f.This critical length is typically of the order of a metre (Knight et al., 1996) at low frequencies.The fine grid cells that we consider are therefore in the range where we may apply Equation 3 while our coarse cells are not.The "patchy" equation (κ f is then the saturation weighted arithmetic average of phase moduli) is recognised as an upper limit for the fluid modulus but it requires that fluids are segregated and may not be fully applicable for our coarse Sengupta and Mavko (1998) have shown that seismic velocities, upscaled from small scale simulations, are closer to the uniform saturation case than patchy.We therefore prefer to use Equation 3 at both scales and calibrate the differences as part of the error.
Finally, the interpolation of predicted impedances down to the fine grid also adds to the model error.Sharp transitions of the fine grid impedances, which may be due to saturation changes at the waterfront or geological variation, cannot be captured accurately.The covariance matrices of the model errors are then obtained from ε i m .We calibrate errors by modelling at two scales similar to O'Sullivan, 2004.If our fine grid model is reasonably accurate (x i t ≈ x i fg , the fine grid result), then we can write the model error as the difference of fine and coarse grid simulation results. (10) The scale change error may consist of a constant systematic term, a parameter related term and a term related to the geostatistical realisation of the geo-model.For this paper we will ignore the parametric contribution, which may occur as different parameter combinations shift the flood front, for example.To first order approximation, we can then use the mean error over a number of realisations to modify the misfit function: (11) In other studies the parametric effect on the error was calibrated (e.g.O'Sullivan, 2004;Glimm et al, 2004) and a similar approach could be applied here if necessary.

Sampling from the Parameter Distribution
For a given ensemble of models, we use the Neighbourhood Algorithm together with existing parameter vectors and misfits to resample the parameter space.Resampling takes place by dividing the parameter space into Voronoi cells.The best n r models are selected and n s /n r new models are randomly chosen in each Voronoi cell in this sub-sample.Large values of n s and n r improve the chances of avoiding local minima and should be increased with the dimension of the parameter space (Sambridge, 2001).

Quantifying the Scale Error
In this section we quantify the errors that arise when we use our coarse scale geo-model as input to the simulation process and comparing it to a fine scale result assumed to be the truth.
We generate a fine grid geo-model (Fig. 3a) by distributing Net:Gross and permeability (horizontally, it is isotropic and independent of Net:Gross).Vertical permeability is a function of Net:Gross and porosity of the sand is uniform, as in our Schiehallion field study (Stephen et al., 2006a).Properties are distributed using sequential Gaussian simulation using GSLIB (Deutsch and Journel, 1992).We generate multiple realisations of the permeability with values at the wells fixed from the first realisation.We fix the Net:Gross assuming that it has been constrained by 3D seismic.transmissibilities should be used in the coarse grid simulation model rather than permeabilities, which are merely estimated in this method.For our models the permeability estimate is very good, however, and we use that here.Simple volume averaging is used to obtain the coarse grid Net:Gross (Fig. 3b).

TABLE 1
Default petro-elastic transform parameters for the dry bulk modulus in Equation 1, obtained from history matching (Stephen et al., 2006c) Following two-phase simulations (Eclipse Manual, 2006), the petro-elastic transform is used with the parameters in Tables 1 and 2 to calculate maps of impedance for each model at both scales.The coarse grid impedances are interpolated so that comparison can be made at the fine scale.Figure 4 shows the impedance change after one year of production for fine and coarse scale models.The impedance is reduced around the injector due to pore pressure increase but increases as water invades the pore space, thereby stiffening the medium.Qualitatively, the match is quite good although some fine grid details are lost in the coarse grid model.Figure 5a shows a cross plot of the impedance changes shown in Figure 4. Most of the impedances lie on the y = x line and the error appears at the edge of the waterfront (Fig. 5b) where the saturation changes sharply.Impedance changes relative to the base line are dominated by the pressure change and the fractional error is relatively small (Fig. 6a).The fractional error increases for differences of monitor surveys, however (Fig. 6b) when the flow rate remains almost constant and so pressures changes are small.The impact of the error at the waterfront is more apparent.The scale error appears to be predominantly a function of the saturations.When we apply a saturation dominated petroelastic transform (E κ = 0 and fluid moduli are independent of pressure) at both scales we obtain a similar result (Fig. 7a) but when we make the stress sensitivity of the dry frame quite large (P κ from MacBeth, 2004), the apparent scale differences are quite small (Fig. 7b).
To calibrate the error, we run a number of simulations with different realisations of the fine grid geo-model.Each geo-model is upscaled and the impedance error between the coarse and fine scale is calculated.Figure 8    Impedance differences (first minus second year) as (a) cross-plot of fine and coarse grid and (b) map of differences.
Figure 7 Cross plots ofiImpedance differences (pre-production minus first year) for the geo-model in Figure 3 where the petro-elastic transform parameters have produced a response that is a) saturation dominated and b) pressure dominated.
first and second years of production.The error is systematic and sizable and reflects the interpolation effect due to sharply changing impedances at the waterfront.The mean error is then included in the misfit calculation (Equation 11).
The contribution to the errors from the simulator can be estimated by upscaling saturations, pressures, porosities, etc., from the fine grid simulation and applying these on the coarse grid when we calculate impedances.The errors that remain will be due to the effect of applying the saturation law (Equation 3) and those from interpolation (Equation 5). Figure 9a shows the average model error when we do this for the first year of production.In Figure 9b, we take the difference between the original error and this reduced mean error to estimate the contribution from the simulation alone.We see that there is a sizable contribution from the simulation error, particularly between the waterfront and the injector well.
With a calibrated error we can then modify the misfit function and begin history matching.

History Matching
In this section, we determine the impact of the model error on history matching examining situations where it is absent, ignored or is compensated for.We also compare it to the data error to determine relative importance.
Figure 10 shows a coarse grid horizontal permeability distribution, which has been obtained by upscaling a fine grid geo-model as in the previous section.The coarse grid geomodel forms our base case in history matching and is modified by changing permeability multipliers at the pilot points shown in Figure 10    Misfit comparison where a coarse grid geo-model has been modified using historal impedances from: base case coarse grid (Case I, black), the fine grid (Case II, red), error corrected fine grid (Case III, green), fine grid plus typical noise (Case IV, yellow) and error corrected fine grid plus data error (Case V, blue).The base case coarse grid geo-model was obtained by upscaling from the fine grid and the solid lines indicate the misfits of the base case for II-IV.pilot points.We populate the initial ensemble of models by randomly sampling a six dimensional parameter space such that the pilot points multipliers can vary from 0.1 to 10.The initial ensemble is independent of any knowledge of the base case geo-model.The fine grid representation of the geomodel is used to generate our "truth" impedances.
We also apply synthetic data errors in Equation 11 based on the statistics of the observed impedances from the Schiehallion field.After suitable conversion, the field study data error is equivalent to a standard deviation of impedance of 0.015 g/cc km/s for each time lapse difference map.Correlation in the field data error was negligible at the coarse scale so we generated noise in the fine scale impedances that was Gaussian and uncorrelated.
To estimate the effect of the model error on seismic history matching we modify the same coarse scale geo-model by matching to truth impedances from: -A coarse scale geo-model directly upscaled from the fine grid realisation.This is a simple way of determining how history matching performs in the absence of errors.The parameter error may be reduced to zero completely in this case (e.g.Stephen et al., 2004).-A fine grid geo-model from which the base case geomodel is derived by direct upscaling.The mean model error is ignored.-As (II) but now the mean model error from Section 2.1 is included in Equation 11. -As (II) but now noise is added to the fine grid observed data to emulate noise observed in the field study (Stephen et al., 2005).We use the data error to derive an uncorrelated random noise field, which is then added to the observed impedance map.-As (IV) except that the mean model error is included in Equation 11.
Figure 11 shows a comparison of the evolution of the misfit for each of the 5 history matching cases above.In Case I, the misfit declines at a log:linear rate with increasing ensemble size and has reduced by two orders of magnitude after 1024 models have been simulated.In Case II, the misfit reaches a plateau due to the model error and we are close to convergence.Note that the misfit of the base case geo-model is greater than the global minimum misfit.When we include the mean error in Equation 11, Case III, we remove a significant and systematic error such that our misfit continues to decline and convergence requires many more models.Once again, however, the base case can be improved upon.In Case IV, the data error affects both convergence and the best model in much the same was as the model error.We can correct the misfit to compensate for the error to a reasonable degree in Case V but a plateau effect remains.In a real field, we would not know if the model error is important without such fine grid simulations.
We analyse the uncertainty of the parameters by converting the misfits into a Posterior Probability Density (PPD) via Bayes Law (e.g.Stephen et al., 2006b).We can then resample the PPD using Markov Chain Monte Carlo Methods (Sambridge, 1999b) to obtain distributions of the permeability multipliers.Figure 12a shows the distribution when there are no errors present.The effect of errors can be seen in Figure 12b such that the peaks are shifted.By including the mean model error in Equation 11, we get a result that is closer to the zero error case (Fig. 12c).We have shown that the model error affects misfit convergence and is equivalent to the data error in size and effect.It is possible to remove the bias, which improves the history matching process and gives a better estimate of the most likely model and its uncertainty.

DISCUSSION
In this paper we have identified the scale dependent errors that arise from the various modelling steps taken to predict changes in seismic impedance.As we coarsen from finer to coarser scale simulations, the waterfront becomes more spread out.The saturations are then used in a petro-elastic transform containing a saturation law that is inappropriate for the scale.Finally, interpolating predictions to the measurement scale creates errors of representivity.
In this study, Figures 7 and 9 indicate that the error is largely located at the flood front but that may not always be the case.The saturation error can be reduced or the pressure error increased as a function of the reservoir parameters.The former occurs if there is a high degree of heterogeneity at the sub-grid scale.In that case, the numerical dispersion error, apparent from Figure 9a, may actually be comparable to the physical dispersion induced by the heterogeneity.Alternatively, the relative permeability curves may have captured the effect of sub-grid heterogeneity with appropriate upscaling.The waterfront would then be spread out naturally as well as numerically and the coarse scale simulation will then be closer to the truth.The saturation error is reduced when we increase the stress sensitivity because the spatial variation of pressure is more gradual.In cases where we have sharp changes in permeability, such as at faults or at the edges of high permeability zones from facies transitions, we may still have a large model error, however.This would be the case, particularly, if permeability transitions are not co-incident with the edges of grid cells.The coarse grid should be defined relative to the faults before the geo-model is built and this is often the case.It is harder to avoid intra-coarse grid transitions of permeability during facies modelling, however.To do so requires that the facies transitions be modelled deterministically, perhaps constrained by 3D or even 4D seismic.
We calibrate the error for a simplistic case where only permeability is modified.The permeabilities are generated using sequential Gaussian simulation with an assumed known range for the variogram.In reality, this may have some uncertainty, bounded by the size of the geo-model and the grid cell size and the geologist's input.We may then choose a single realisation as a base case to be modified.It is extremely unlikely that our modified geo-model will match reality and so several realisations should be considered to reduce the inherent parameteric error.In addition, the Net:Gross variable affects the balance between saturation and pressure effects in the petro-elastic transform such that alternative realisations could change the model error.We can constrain Net:Gross to an areal map in many cases but the fine scale is often uncertain and an appropriate modelling scheme is required.Again the realisation may be important.Finally, the error may depend on the parameters of the system that we vary during history matching.If we move a permeability boundary or force the waterfront to a different location, the model error will change.O'Sullivan ( 2004) produced a method of modelling and parameterising the errors and interpolating as a function of the unknowns.A similar scheme may be necessary in general.
The error in our models has significant impact on the history matching result and is unavoidable without seriously increasing the detail that we must consider, which in turn would reduce our ability to search the parameter space sufficiently.We can accept the error provided that we that we compensate for it in the misfit calculation so that we have a better measure of how close our predictions are to reality.We show that it is possible to correct for the errors in the modelling shown here so that we find more accurate models.Of course, the effect of the model error must be considered relative to the data error, which, if large, may dominate the history matching process.By including noise of similar magnitude to that observed in the Schiehallion field, we find that we reduce our ability to find better models.If the level of data error is much larger than we have observed, it may not be worthwhile spending time calibrating the model error.

CONCLUSIONS
• Scale dependent simulation and petro-elastic transformation errors have been identified in seismic history matching and are of equivalent size.• The model error produces a non-zero minimum misfit and affects the accuracy of history matching.• We can calibrate the error and reduce its effect on history matching.• The errors dominate at sharp changes in the impedance field, either at the saturation front or at faults but saturation errors are apparent elsewhere.• If faults or other permeability changes are absent, a pressure dominated seismic response will have negligible model error.
Figure 1Schematic of the iterative automatic history matching process where many models are generated per iteration and their misfits are used to pick new parameter values.
Figure 4Maps of calculated impedance differences (baseline minus survey) for the geo-model in Figure3after the first year of production for a) fine grid geo-model and b) upscaled equivalent after interpolation.

Fine
Figure 5 Impedance differences (pre-production minus first year) from Figure 4 as (a) cross-plot of fine and coarse grid and (b) map of differences.

Figure 6
Figure 6 Figure 8Average model error of impedance change calculated for six realisations of Net:Gross, porosity and permeability constrained to the same well data as in Figure3for differences (a) pre-production first year and (b) first year minus second year.

Figure 9 a
Figure 9 a) As Figure 8a except that fine scale saturations, pressures etc. were upscaled to the coarse grid and then the impedances were calculated; b) Figures 9a subtracted from 8a to estimate the simulation error.

Figure 10 Example
Figure 10Example permeability model for the coarse grid with pilot points shown.

Figure 11
Figure 11 Figure 12One dimensional marginals of the pilot point multipliers used to modify the permeability distribution for a) Case I where the coarse grid impedances are matched; b) Case II where the fine grid impedances are matched and c) Case III where the average model error is used in Equation11.The numbers in the legend correspond to the points identified in Figure10.The base case multipliers are unity in the base case (i.e.truth).
Batzle and Wang (1992)transform parameters for Gassman's Equation, 2 and the saturation Equation3.The fluid bulk moduli are functions of pressure, temperature etc. and full details of the calculation can be found inSoldo (2005), which followsBatzle and Wang (1992)