The pre-drill and real time while drilling overpressure prediction based on effective basin model concept Текст научной статьи по специальности «Науки о Земле и смежные экологические науки»

The pre-drill and real time while drilling overpressure prediction based on Effective Basin Model concept

A.G. Madatov, A.-V.I. Sereda

Higher Mathematics Chair of the Polytechnical Faculty, MSTU

Abstract. The problem of the pre-drill predictability of an overpressure phenomenon has been considered within the range of the Effective Basin Model concept. The inherent and numerical origins of the relevant prediction scenarios have been analysed. The solutions of numerical problems are attached in the form of the appendixes. The approach to the stabilisation of an inherent prediction error has been proposed.

The real time correction of a pre-drill prediction in during drilling ("Look Ahead") mode has been considered in the context of the re-calibration or/and adjustment of a prediction model and a framework earth model respectively. The re-calibration problem has been introduced as an extension of the single-well data inversion problem. The corresponding ill-posed task has been formulated and the approach to getting its locally unique solution has been discussed. The relevant computer algorithm has been described in details. The applicability of the new approach and solution to the real time drilling conditions have been tested based on the speculative (North Sea area) and real data (Gulf of Mexico) examples. These tests confirm high converging speed of the updating routine and flexibility of the algorithm. The achieved results correspond to higher accuracy in comparison with available standard methods.

1. Introduction

Pore pressure estimation and prediction is a difficult and challenging data analysis problem. More so this is true for geologically old sections, which have had passed through a complex sedimentation history (Swabrick et al, 1999). It is vitally important to get an earth model with an adequate level of complexity to establish a prediction framework relevant to data available a priori and during drilling (Jackson, 1979). From this point of view a combination of basin modelling approach with data inversion technique introduced in our previous papers (Madatov et al., 1997; 1998; Madatov, Sereda, 2000a) looks as the most promising strategy in connection to the overpressure prediction problem. This viewing coincides with the general conclusion recently coming from industry practice (Duppenbecker et al., 2002; Nordgard Bolas et al., 2004; Williams, Madatov, 2005; Standifird, Matthews, 2005).

The Effective Basin Modelling approach allows forward modelling of pore pressure in basin time scale with all practically important 3D phenomena accounted (Madatov, Sereda, 2005a). It is important to stress that a rough estimation of the overpressure effects within a frontier area can be performed prior drilling just by using the fundamental link between retaining and generation geo-fluid potentials developing in sedimentary rocks during their burial history and associated processes: compaction, heating, diagenetic transformations and so on. Still such a model cannot pretend on the pre-drill prediction within a virgin region, with no offset wells associated in tectonically common block with target well location. It only allows extension of data driven seismic methods (Dutta, 2002) by adding geological information normally available a priori such as age, thickness and lithology of the sediments. The situation changes there and then where and when some offset well data can be involved for EBM calibration1. In this case a multi-well data inversion in conjunction with prior geology knowledge allows proper calibration of the framework earth model established for the calibration area where offset wells exist (Madatov, Sereda, 2000b). Pore pressure can then be predicted away from the calibration locations based on variations in geology plus interpolation/extrapolation of control model parameters.

Evidently, this scheme describes an approach to a newly developed strategy of overpressure prediction very roughly. The series of articles published by us in the MSTU Proceedings earlier was dedicated to development of an fit-for-purpose basin model (Madatov, Sereda, 2000a; 2001; 2005a); the model upscaling (Madatov, Sereda, 2003) and calibration of the model based on a single well (Madatov, Sereda, 2000a,b) and multi-well data inversion (Madatov, Sereda, 2005b). The given article is aimed to summarise the achieved results in getting practically applicable strategy of overpressure prediction prior drilling with possibility for timely updating it in real time during drilling.

Below in the text we will use axiomatic and terminology introduced by us in above mentioned papers. Some additional general level terminology related to the math and the prediction topics is available in Appendix 1.

1 It is important to pick out that the overpressure problems are normally associated with frontier sub-blocks and/or deeper horizons within the well explored areas. The exploration drilling hazard in the Gulf of Mexico deep water is a typical example of these circumstances.

2. The predictability of the overpressure phenomenon and generic strategy of its prediction based on the EBM approach

The predictability is important feature of an Earth model grounded on Model Recognition principle, which has been introduced by Braverman (1962) as the following: "The number of sensitive to the target phenomena properties of a proper prediction model is finite and relevant set of sensitive model parameters can be recovered based on model calibration".

In a very general manner the predictability of any geo-phenomenon including overpressure implies satisfying to three following conditions (Slingerland, 1990):

1. There is at least one earth model capable to explain and to describe quantitatively a target phenomenon based on a priori knowledge, available empirical-deterministic basis and intuition.

2. This model is recognisable based on available data.

3. The Prediction Model can be established, i.e. the sensitive descriptors can be continuously extended (projected) into the target location in space and time coordinates from reference locations within the calibrated earth model.

As it follows from this definition the key factors controlling predictability of a natural geo-phenomenon are suitable for calibration earth model attributed with the minimally sufficient number of stably consistent and continuously sensitive to target model parameters. The concept of effective basin model introduced in (Madatov, Sereda, 2005a) gave the required earth model in the form of rather general combine 3D hydrodynamic (HD) basin scale model, which was attributed by a finite number of the sensitive to overpressure model parameters, which could be reconstructed based on inversion of pressure related data.

Computer engine for well interpolation and pressure modeling

0

Target well

■;.75 2 2.25 Pressure, [g/cc]

Fig. 1. The scheme of pre-drill overpressure prediction based on prior calibrated EBM model

a - components of the Prediction Model in the map form; b - the Prediction Model in the form of the stack of the calibrated model parameter maps tied to the offset well positions (marked with the triangles); c - computer engine for well interpolation and pressure modeling; d - overpressure prediction scenarios for target well position

It is important to pick out that most of the key overpressure mechanisms are reproducible within the range of this modelling basis. The model descriptors are given here in the form of lithology depended on control model parameters assigned for every formation element according to the model parameter space specification (see Madatov, Sereda, 2005a). The multi-well data inversion routine provides the means for the calibration of the established within the area framework earth model (Madatov, Sereda, 2005b). Thus, the predictability of overpressure can be ensured based on the EBM concept. The quality of overpressure prediction here can be quantitatively a priori and posterior controlled based on analysis of every local solution of the data Inverse problem (its uniqueness, stability and consistency with available geology knowledge).

The_genfiric„Siraifigy.flLflVfir:pressure„predictjan follows from the above introduced predictability definition. Regardless of the basis earth model and set of sensitive model parameters it implies consequent fulfilling of the following sub-problems:

• Getting of a Framework HD Model to include: setting of a background HD basin model suitable for accommodation of the prior knowledge about local geology, fluid dynamics and its upscaling against measurable in situ geo-fluid system properties and the relevant geophysical data.

• Getting of a Prediction Model to include: calibration of the tuneable model parameters attached to a framework HD model against definitive multi-well data set by using specific inverse problem (SIP) solver in the form of the e-equivalent parameterisation sets.

• Getting of Pre-drill Pore pressure scenarios to include: interpolation of the required model parameters at the target location within the framework model for all e-equivalent prediction models and computing synthetic pore pressure curves by using the relevant forward problem solver.

• Getting of a Definitive Prediction curve to include: ranking of pore pressure models based on stability and consistency analysis. Choosing the multi-curve average result for the highest probability rank of the prediction model as the most likely scenario for pre-drill prediction.

The sketch in Fig. 1 illustrates the main phases of generic prediction strategy: calibration (see Fig. 1a); build up a prediction model (see Fig. 1b); producing pre-drill prediction scenarios (see Fig. 1d).

The sub-problems 1 and 2 of this list were matter of some previous articles (Madatov, Sereda, 2000b; 2001; 2003; 2005a; and Madatov, Sereda, 2000a; 2005b, respectively). A brief summary about them could be formulated as the following.

Setting and upscaling of a framework effective basin model

The combined approach to the 3D forward problem for HD system in the basin time scale allows construction of the numerical solution in the form of consequent combining 1D and 2D problem solutions at the common time grid. Thus, the relevant model specification is representable in the form of the stack of the 1D (available offset well models) grids tied to the stack of 2D grids (subsurface geometry for available aquifer formations).

Setting of a framework earth model includes: establishment of common for the calibration area stratigraphy column at formation level of integrity; assigning empirical relationships for porosity-permeability relationships for every litho-type available; successive specification of every formation unit in agreement with model space specification rules.

Upscaling of an earth model available at the previous step implies optimal simplification of the formation content and fault geometry model.

Framework model calibration based on developed inverse model solution approach

The depth-lithology depended real data quality estimation is an inherent feature of the given framework

model.

The final accuracy of calibration data causes a set of equally valid prediction models to be formally fitted to the available data sets. In other words, it is inherent property of non-unique IP solutions in practice.

The SIP solution based on multi-well data set for the framework model provides formally unique calibration result built on posterior satisfaction to the additional to data quality RV and Sim criteria.

The more a priori restricted is domain of potentially suitable sub-class of HD models within the model parameter space the more unique calibration and more certain prediction are.

Speaking about inherent non-uniqueness of the calibration results it is important to distinguish between local and global non-uniqueness of the relevant multi-well data inversion problem. Formally all achievable SIP solutions deliver satisfactory misfit between real and synthetic pressure-porosity data at any and all of the offset well positions in the form of e-equivalent prediction models. Still, the locally non-unique solutions form simply-connected compact sub-domain Xj within the model parameter space X which can be interpreted as a location of some single sub-class of HD models in there. Contrary to it, globally non-unique solutions belong to different local mínimums interpretable as compact sub-domains {XJ} having empty crossing space.

Scanning through all of the priory predefined sub-class of HD models {Xj} allows generation of a certain sequence of SIP solutions {xj e XJ}, which are all e-equivalent in terms of satisfying to the common data stop criterion and generally non equal in terms of satisfying to posterior similarity (Sim) criterion and variability (RV) criterion (Madatov, Sereda, 2005b). Every particular solution achieved for the established framework model in the form of its calibration provides establishment of the relevant particular prediction model for the calibration area. The quality of every prediction model could be then non-formally estimated basing on prior expert analysis of the predefined for calibration sub-class of HD models and also basing on posterior interpretation of the model parameters' distribution for the particular prediction model.

Evidently, the terms interpretation and expert analysis carry implicit connotations of uncertainty. An ability to quantify this uncertainty in the context of the overpressure phenomenon's predictability provides logical basis for drilling risk assessment. Usage of the Effective Basin Model concept (Madatov, Sereda, 2005a) and formal criteria elaborated for SIP solutions (Madatov, Sereda, 2005b) obviously protects us from creation some geo-science fictions, which have nothing in common with real nature. Still, the quality of produced predictions in general and the definitive overpressure scenario in particular requires more rigorous analysis of possible errors to be properly assessed (Hunter, Mann, 1992).

Beside of such obvious steps at any prediction technique like cross-checking of results wherever possible, blind tests etc. there are at least three specific directions for minimising drilling risk to be addressed in connection to the "pore pressure surprise" (Brekke, Thompson, 1996).

The first direction is non-formal evaluation of consistency between any of the prediction models and available a priori results of particular sedimentary basin analysis (Lerche, 1990; Miall, 2000) and/or the relevant deposystem forward modelling (Merriam, Davis, 2001). Consideration of this approach is a matter of separate study, which results are far beyond the scope of the given article.

The second direction is to emphasize, and to quantify wherever possible, the uncertainty and stability associated with every prediction model in order to pick out the most stable and certain prediction as a definitive overpressure scenario.

The third direction is to make sure that the pre-drill prediction is converging to post-drill overpressure estimation while real time prediction updating, i.e. in during drilling mode.

Last two directions are illuminating in the given article. In particular, in the below sub-section the origins of pre-drill overpressure prediction errors are analysed in the context of the introduced above generic prediction strategy.

2.1. The formal definition of pre-drill overpressure prediction scenarios

The pre-drill overpressure prediction scenarios wherever inside of the calibrated Framework Model are getting by multiple runs of the forward pressure modelling for all available in calibration area Prediction Models.

Each particular prediction scenario can be represented in the operator form as the following:

p + Ap = M [(z + Az); x + Ax], (1)

where the vector p together with its Ap vicinity are defined on the data comparison space P in agreement with Data Space specification (Madatov, Sereda, 2005b) in the form of the most likely curve with the error bars attached (see Fig. 2e-f); the vector x together with its Ax vicinity are defined within a priori bounded compact sub-space on the model parameter space XJ c X as a set of e-equivalent solutions of the inverse problem for calibration multi-well data set (Madatov, Sereda, 2005b); M[ ] is the final difference operator implemented for solution of the 3D HD forward problem in the basin scale defined on the model parameter space X together with the boundary conditions (Madatov, Sereda, 2005a), Note that the initial conditions of the relevant forward problem is getting from the available sub-surface geometry-stratigraphy features established for the relevant framework model. They are represented explicitly in the general formula (1) in the form of the vector z together with its Az vicinity defined in agreement with the framework setting and interpretable as a possible variation in the formation depth and/or age within the calibration area (the compact sub-space Z).

2.2. The nature and formal definition of pre-drill prediction errors

Let some j-th sub-class of HD models xj e XJ with an expert-based credibility rank Rj attached have led to successful solution of the locally unique calibration problem and the set of the relevant prediction models now can be established within the area. Let also the target position for prediction be established inside of the calibration area. Let also the framework model be free from error, i.e. Az = 0 in formula (1). As it follows from the theory of multi-well data inversion (Madatov, Sereda, 2005b) all prediction models generated from calibration results are equivalent in terms of model-data misfit accuracy, or just are e-equivalent. The total

number of such e-equivalent prediction models produced from the same start point xj0 inside of the given model type Xj does not exceed the number of the offset wells. Assuming that there is more than just one single option for start point to be assigned inside of the given model sub-class XJ; one can state that amount of e-equivalent prediction models K could be reasonably big for producing statistically consistent estimations.

Thus, the most likely overpressure scenario ea for the target position inside of the calibration area can be produced by means of multi-curve weighted averaging over all synthetic pressure curves computable in agreement with generic prediction strategy. Then the relative pre-drill prediction scenario error is representable in the following form:

|| Ap || /1| p" || = E-[ ||M(z, Xjk) - M(z, Xj") || / ||M(z, Xj") || ] + ¡M, VXjk e Xj, (2)

k = IK

where E{ } is the operator of weighted averaging; M[ ] - the final difference operator of the HD forward problem solution defined on the model parameter space X together with its inherent numerical stability estimation ¡M; the vector xj" represents the most likely control parameter set at the target position given by averaging over K available e-equivalent prediction models produced for the given model type XJ; the vector xkj represents current k-th control parameter set at the target position achieved for the particular j-th prediction model.

In the general case the number of sub-class of HD models can be finite and fixed as J. The relevant generalisation of the pre-drill EBM based prediction errors is defined as follows:

||Ap||/|| fi"|| = E{ Rj [E- ||M(z, Xjk) - M(z, XjA)||/||M(z, Xj")||] + ¡M, VXjk e Xj, (2*)

j = 1,J k = 1,K

where Rj - the weight assigned for j-th type of prediction model in agreement with its credibility rank and the rest notifications are defined above.

Evidently every particular model parameter vector xj defined at the target well position can be represented in the form:

xk = Xj" + Axkj. (3)

The subset {Axk}j defined over all (k = 1,K) particular prediction models established for the given model sub-class XJ can be understood as a result of combining of interpolation errors Zix with uncertainties 8kx of the current e-equivalent prediction model. Here lower indexes i and k indicate the target point location within the calibration area and the particular calibration result achieved for the given model sub-class respectively.

Let the position of a target well inside of the calibration area be fixed. Thus, the interpolation error is certain as Zx. Then local vicinity of most likely vector xj" position can be estimated in average as follows:

Ax = E-{8xjk + Zxjk}, Vxjk e XJ. (4)

k = 1,K

Note, that generalization of formula (4) for all acceptable during calibration model sub-classes {XJ} implies aggregation of elements {8xjk + Zx/} in the global weighted averaged sum along j = 1,2,. .J given by analogy with (2*) in the form:

Ax = E- {Rj [E-(8xjk + Zxjk)]}, Vxjk e Xj ; VXj c {Xj}. (4*)

j = 1,J k = 1,K

Every k-th solution xjk of multi-well data inversion problem delivers a conditionally unique e-equivalent solution for every involved offset wells (Madatov, Sereda, 2005b). It means that the data fit quality does not exceed certain misfit accuracy e predefined for all offset wells. Since, SIP solution is stabilized (see Madatov, Sereda, 2005a for more details) one can claim that this statement remains true wherever inside of calibration area on condition that the framework model or/and triangle mesh do not generate numerical instability in forward modelling or/and interpolation routines respectively.

Note, that except numerical instability of an interpolation routine the interpolation error is several orders lower than the data»model fit error accepted as the e-equivalent SIP solution (Zxkj << 8xkj Vj = 1,2,..J). Thus, the sub-index "k" at the term Zxkj in formulas (4-4*) might be omitted. Yet, we have also to include formally in consideration instability ¡M of the final difference forward modelling operator M[ ] as a part of inherent numerical modelling instability. We have explored both interpolation and numerical unstable situations and results are presented in Appendixes 2 and 3 respectively.

Note finally that on conditions of the zero tolerance error (Az = 0) introduced before for the framework forward model and beside of the mentioned unstable situations the relative pre-drill prediction error (2-2*) does never exceed misfit accuracy e predefined as an area common stop criterion.

What will happen if Az ^ 0 in formula (1)?

2000

3000

4000

5000

1,0 1,2 1,4 1,S 1,8 2,0 2,2 Overpressure [g>c£]

a

1000

2000

3000

4000

5000

1 ,0 1 ,2 1 ,4 1 ,6 1 ,8 2,0 2,2 Overpressure [g/cc]

b

I

1 1 v

J

\

M

1 □ 1 2 1 4 1 6 1 8 2 □ 2 .2

Overp ressure

e

[g/cc]

l

■ \ ; p

I i.

- - c a C \

1 u-L.

P

1,0 1,2 1,4 1,B 1,0 2,0 2,2 Overpressure [g/cc]

f

Fig. 2. Evaluation of the pre-drill prediction quality based on two HD model sub-classes used at the

prediction model calibration

a,b - the pre-drill prediction scenarios evaluated for the target well position (marked on below maps with the star) from the set of prediction models calibrated against subclass 1 (a - active HC generation from source rock) and subclass 2 (b -mud-rock permeability drop). The code of curve corresponds to base well number at the particular SIP solution; c,d - the pre-drill prediction in the form of pore pressure map tied to the target formation "F3" (highlighted on plots a,b,e,f); e,f -the definitive pre-drill prediction scenarios for prediction model 1(e) and 2(f) with error bars attached

d

c

The example in Fig. 2 illustrates the situation where two pre-drill predictions were made basing on two different prediction models. The first prediction model x (Fig. 2a,c,e) was calibrated for the EBM sub-class Xi with the active fluid generation from source rock predefined as a key overpressure mechanism (Xiaorong, Vasseur, 1996; Gandi, Berg, 1997) (see also Fig. 3a). The second prediction model x1 (Fig. 2b,d,f) was calibrated for the EBM sub-class X2 with mud rock permeability drop predefined as a key overpressure mechanism (Magara, 1978; Osborne, Swabrick, 1997) (see also Fig. 3b). The definitive pre-drill scenario curves were produced basing on the weighted multi-curve averaging according to routine (2), where every particular prediction curve was computed for the same target position as follows

pjk = M [(z; xjk + Axjk], j = 1,2; k = 1,2,.. K (5)

The multi-well data set have included 12 offset wells, i.e. K = 12 (see also Fig. 5).

Two formal confidence intervals built based on 12-curve summarizing appear to be rather different (compare Fig. 2a,e with Fig. 2b,f respectively). The narrower confidence interval is related to the second prediction model (X2), whereas some of the prediction p1k based on the first prediction model (X1) doesn't fit to the definitive prediction p1" with the acceptable error despite of the fact that all of the included in multi-curve averaging prediction models deliver e-equivalent SIP solutions for all offset wells.

We suggest the term Apparent Modelling Error at the target position for recognition of this prediction instability case. The explanation of the relevant numerical feature lies in the insight into the inherent stability of the given sub-class of HD model to the small variation Az of the framework model vector z defined for the target position. This kind of modelling error apparently appears due to instability of the given model type against the framework model variation and tolerance errors in it.

2.3. Apparent modelling errors at the target position

Let the framework model be fixed in terms of the vector z components (see description at formula (1)) and also be calibrated based on the given model sub-class Xj. That means in particular that the components of the model vector x are defined at every offset well position to be delivering stable and conditionally correct e-equivalent solutions of SIP (Madatov, Sereda, 2000a). The rest of components required according to the model specification for the whole area including the target location are getting recovered by interpolation between calibrated components (see Fig. 1b,c). This process is routinely repeating for K start points of multi-well data inversions (Madatov, Sereda, 2005b). Finally K prediction models {x} are formed within the area and consequently K prediction scenarios can be computed for the target location in agreement with (1). Note that every particular prediction model satisfies to the data stop criterion at every offset well, since e-equivalent SIP solution was reached at every particular calibration case.

The question rises whether the sub-set (3) recovered by interpolation on any target position inside of the area and then mapped on Data Comparison space P by using the continuous operator M [z;{Ax"}] will also form compact sub-set {ep"} (e-locality) around of the most likely data prototype p" (prediction) given by multi-curve weighted averaging routine (2)?

The vector solution x of the multi-well SIP satisfies to the minimum variability criterion. This provides minimal area variability ranges for all involved in tuning control model parameters (non-fixed components of the model vector x) as it was proven in (Madatov, Sereda, 2005b). Consequently the recovered on the target well position components of the vector x also will not exceed the constrained variability ranges on condition of the stable interpolation routine. Thus, control model parameters on the target location will belong to the relevant A-vicinity of the vector solution x" as soon as all particular e-equivalent solutions are locally unique. As the modelling operator M[z; x] is numerically stable on the model space X where it is defined and the compact sub-set {Ax"}c Xj is simply-connected, it can be uniquely mapped on Data Comparison space P in the form of the corresponding subset {ep"}. Thus, the prediction quality seems to be constrained by the predefined data accuracy criterion.

However, the sub-surface geometry at the target position generally differs from the ones at the offset wells. This creates prerequisites for apparent modelling errors, which could corrupt definitive prediction quality. The cause of this is a sub-surface geometry variation and/or tolerance error potentially available in the framework model. Indeed, the stability of mapping result M[z;{Ax"}] is now depended on how stable is the forward modelling operator M implemented wherever inside of the given model sub-class Xj to the area variations in the framework model vector z. We address readers to the Appendix 4 where the local numerical features of the combine 3D forward modelling operator have been analysed. Here we just stress that it is highly non-linear. Thus, the small variation in each or in both of input vectors z and x can significantly disturb result of mapping on Data Comparison space P.

In the context of the considered calibration problem it was shown that the prediction quality could be ensured wherever inside the calibration area on condition that the vector z remains non-changeable, which is often not true. In more general and practically more important case the sub-surface geometry (components of the vector z) is changeable from point to point within the area, and in particular at the target position it differs from the offset well geometry.

To be more specific we have to consider the calibrated model sensitivity response on perturbation of vector z. We address readers to our article (Madatov, Sereda, 2003) to get more details about sensitivity analysis and sensitivity features of the implemented forward modelling operator. Here we just will determine the sensitivity vector and will give an important conclusion about stability of mapping result M[z; x] without letting it mathematically proven.

The sensitivity vector is defined at any arbitrary position of the model parameter vector x e X and the framework model vector z e Z as follows:

S(x,z; dx,dz) = ||M[z + dz;x + dx] - M[z;x] || / ||M[z;x] ||, (6)

where dx,dz are the infinitesimal local perturbations of the vectors x and z such that || dx ||" || x || and ||dz||" ||z||; the forward modelling operator M[z;x] = p is defined above; ||M[z;x] || is the L2 norm of the synthetic data vector p e P is given in agreement with the data vector space P dimension.

Now the following statement can be formulated.

The forward modelling operator M[z;x] provides absolutely stable mapping wherever inside of the local sub-sets {Ax} c Xj and {Az} where the vectors z; x are varying on condition that the relevant sensitivity vector (6) remains monotonously continuous wherever within the sub-sets {Ax} and {Az}.

If the components of the input vectors xA and zA at the target position are getting by interpolation between calibration wells the instability of prediction, i.e. excess over the data stop criterion can only be result of two possible prediction cases:

Prediction cases 1: The prediction model at the target position is outside of the given sub-class XJ. Prediction cases 2: The monotonous continuity condition for sensitivity vector is not satisfied inside compact sub-spaces XJ and Z.

The first case may correspond to the situation when the prediction target belongs to a different tectonic block separated from the one contained all offset wells by a tectonic fault. It is clear that such case is recognisable in advance based, for example, on the structure interpretation of the seismic data. In the general case, based on the interpretation of the calibration results it is always possible to subdivide the framework area onto sub-areas or tectonic sub-blocks inside which the conditional uniqueness of the SIP solution is ensured. One of the possible strategies can be based on getting an area representative sensitivity image of certain EBM subclass (Madatov, Sereda, 2003).

Since the model sub-class Xj is assumed to be fixed we focus on the second prediction case in this discussion.

This case is less predicable in advance since the target can belong here to the same sub-class of the model, but pre-drill prediction could still suffer from instable mapping. The introduced above statement allows controlling accuracy of pre-drill prediction quality if this apparent modelling error occurs basing on prior sensitivity analysis of a prediction model.

This can be illustrated by the example with two model sub-classes introduced above (see Fig. 2).

We have fixed two prediction models x1 e X1 - source rock generation controlled (see Fig. 3 a) and x2 e X2 - mud-rock permeability controlled (see Fig. 3b) for the same target position z e Z and we investigated the L2 norm variation of the relevant sensitivity vector (6) against perturbation of the vector z. Namely, the definitive framework vector z was consequently being perturbed in order to scan through the discrete list of the framework model errors {0,|z,2|z,3|z,...,8zMAX}, where 8zMAX - the maximal tolerance error, - the sampling step. These errors can be associated with possible mistakes in the top formation depth or thickness or both.

The sensitivity response in the stable case (prediction models x2) is monotonous with the close to a constant first derivative d||S||/d||z||, whereas the sensitivity response for prediction models x1 is non-monotonous with the local maximum around depth error -120 m (see Fig. 3c). The same dynamic keeps repeated through full ranges of x1,x2 variations for K prediction models.

The explanation of the observed apparent modelling errors and related to them sensitivity features for type 1 of the prediction model is rather trivial. According to the kinetic model of source rock degradation - HC generation - this process non-linearly depends on heating rate of the source rock during its burial history, and in particular during its passing through HC generation temperature window (Tissot, Welte, 1978). This rate is in turn controlled by subsidence rate of upper formations and the relevant paleo-thermal gradients. The peak of the

gas generation corresponds to a maximal pore volume expansion and to a minimal HD conductivity of source rock (Xiaorong, Vasseur, 1996). Hence it timely coincide with overpressure occurrence in the corresponding depth interval. Clearly, that the maximum sensitivity response in Fig. 3c for the prediction models x1 indicates approaching to the peck generation conditions. Clearly also that any small errors in formation depth and/or age being combined for upper formations (tolerance error can lead to the not-tolerance prediction error if the model vector at the target position is predicted in close locality to the maximum sensitivity response.

Fortunately the sensitivity analysis of a prediction model can be performed prior drilling.

The described example illustrates approach to the prior drilling test of the apparent modelling error during production of the definitive pore pressure scenario. Namely, the particular modelling errors {Sp} caused by only particular tolerance errors {Sz} in the framework model vector components (z) can be estimated according to the following algorithm:

1. The prediction model at the target well position should be computed as pA = M[zA;xA] basing on most likely vectors for both model parameters (xA) and framework model (zA). The relevant definitive pore pressure scenario is considered to be a most likely case among all available pre-drillpore pressure scenarios {p}.

2. The definitive framework vector z then should be consequently perturbed in agreement with the list of the possible framework model errors {0,|z,2|z,3|z,...,SzMAX}, where SzMAX - the maximal tolerance error, -the sampling step.

3. The updated pre-drill pore pressure scenarios should be computed for the perturbed framework model as the following: pU = M[(zA+Sz);xA].

4. The relative modelling errors |pA - pU|/ pA must be checked against tolerance prediction stability limit.

Fm_7 Fm_6 Fm_5 Fm_4 Fm_3 Fm_2 Fm 1

- HC generation constant HD conduction constant

Fm_7 Fm_6 Fm_5 Fm_4 Fm_3 Fm_2 Fm 1

- HC generation constant

- HD conduction constant

Low

Moderate

High

Low

Moderate

High

> 0.6

e s d e

-200 -1 00

0

100

200

Fig. 3. The sensitivity responses to the tolerant framework model errors

a,b - the inherent distribution of the sensitivity vector components for the HD model sub-class 1 and 2 respectively; c - the sensitivity anomaly indicating possible instability in the pre-drill overpressure scenario based on sub-class 1 of the relevant prediction model

Targetformationdepth error [m]

c

2.4. Discussion

The main advantage of the EBM based on the pore pressure prediction in comparison with the traditional seismically driven approach (Dutta, 2002) is possibility to explore different pre-drill prediction scenarios. Quantitative evaluation of prediction quality is especially important for the potentially non-stable situations.

As it follows from above the introduced prediction strategy not single, but series of prediction scenarios can be prepared and analysed for the same target location within a calibration area. Based on the formal and nonformal prior analysis of alternative prediction models it appears to be possible to range them out in agreement with the relevant quantitative and qualitative ranks. The criteria of such ranking are a matter of extended

discussion exceeding the scope of the given paper. In this context it is important to stress that the finite set of pre-drill scenarios accompanied with specifically defined envelop of uncertainties could be prepared for the target position in advance before the drilling operation will be started. Reservation of the full set of the pre-drill pore pressure scenarios ready for substitution of the current definitive scenario ensures drilling engineers from unpleasant drilling surprise.

The fast and stable real time updating of a definitive pre-drill prediction curve is the next important potential advantage of the EBM based approach to the prediction. The relevant "Look Ahead" strategy can be based on real time verification of the Prediction Model consistency with observed in drilling well geology as well as on constant examination of the misfit between available above drill-bit pressure related data and their synthetic prototype and re-calibration (when necessary) of the Prediction Model. This approach is considered in the next sub-section in more details.

3. A real time updating of the pre-drill prediction in while drilling mode

The Real Time (RT) updating of a definitive pre-drill prediction curve can be considered as a final and a crucial stage of the introduced above generic prediction strategy.

Actually, none prediction can avoid the real time adjustment stage, since pre-drill scenario inevitably contents inherent and apparent errors, which must somehow be corrected as new information about framework model and/or sensitive to overpressure phenomenon are getting available during drilling. In practice, however, the converging real time adjustment of the pre-drill prediction is not achievable within the range of industry standard methods. The reasons for that can be emphasised as the following.

Seismically driven pre-drill prediction:

1. The interval velocities available from the surface seismic data analysis have inherent restrictions in both resolvable and noticeable overpressure derived response. These limitations grow up with a target depth.

2. The more purposely oriented technology of the surface seismic analysis is applied the more expensive and time consuming processing should be.

3. The VSP checkshot acquisition required to get data for constraining the original acquisition velocities turns out even more costly and time expensive than seismic reprocessing.

The pre-drill overpressure prediction problems based on seismically driven overpressure response have been considered in the series of the recently published papers (Dutta, 2002; Huffman, 2002; Madatov, 2005) in more details. Here it is important to note that we do not include seismic time-laps measurements (4D seismic monitoring) to the list of methods because the field of this application is rather different from pre-drill overpressure prediction and RT updating. Indeed, this method is practically applicable to history matching problem formulated for the production time scale within the range of pure 3D field scale reservoir model (see, for example, Brevik, 1999; Landro, 2001; Eide et al, 2002; Vasco et al, 2004). The production time-laps included into dynamics of such kind of the model significantly reduce the dimension of the relevant inverse problem, which is now formulating not for absolute data values, but for their time increments. In other words, the relevant 4D framework models are RT updateable by definition.

Basin model driven pre-drill prediction:

1. The architecture of basin models is commonly grid-based and therefore requires mapping of multiple surfaces and intervals over the entire calibration area. On the other hand, the pressure and pressure related data required for calibration are normally presented irregularly along offset well paths. The production of the necessary data maps is not a trivial part of the calibration routine.

2. The complexity of the commonly applied 3D fluid dynamic forward modelling problem that is solved on a geologic time scale ensures that most of the basin modelling approaches are not suitable as a basis for inversion. Indeed, the 3D solution of the relevant inverse problem on a regular grid requires equally regular data distributions of the equal quality in order to start the calibration. In practice, however, real well data are only available at a few locations in such a grid. In addition, simultaneous optimal adjustment of hundreds of model parameters and restoration of the 3D formation surface paleo-geometry or the 3D paleo-stress field distribution within an area can hardly be imagined in practice. From this point of view, the non-uniqueness of the inverse problem solution will increase dramatically as restoration back through burial history continues. Additionally the proper choice of the currently sensitive to the target attributes of a basin model should be based on the multiple recalculation of the sensitivity vector, which is time consuming operation. Otherwise it becomes very depended on the modeler's vision of the HD system architecture and development as well as on his good guess in rock property definition (Lerche, 1990; Ungerer et al., 1990).

3. The final 3D pressure models often take a considerable amount of time to construct and run. The run times on the order of hours to days are not uncommon. This point reduces real-time updates of the basin model derived prediction during drilling to hardly possible operation (Doligez et al., 1985).

The developed EBM concept and the relevant approach to the prediction can be formally attached to the last group of the methods. However, the special combine design of the 3D HD model (Madatov, Sereda, 2005a) and the developed upscaling routine (Madatov, Sereda, 2003) allows to consider it as an especial modification of basin modelling tolls fitted to pre-drill and real-time prediction purposes. Indeed, the following features of the developed method allow to come to this conclusion:

1. The framework HD model doesn't require too much time and expertise from a user to be established and calibrated. (It is counting on a drilling engineer, not a geo-scientist having about two weeks training course as a special background);

2. It is simple but comprehensive enough to reproduce all major overpressure controlling factors in synthetic porosity and pore pressure profiles;

3. The pre-drill prediction model provides accommodation of a multi-source data input normally available during drilling (formation top/thickness estimations, wire logs, overpressure shows, drilling tendencies and so on);

4. It is updateable by inversion against accessible well data in practically reasonable time (during minutes by using laptop-like equipment)

5. It gives possibility to check all available prediction scenarios and getting current definitive decision for both a beginner (by using fully auto mode approach) and advanced experts (by delegating him control on ranking of each scenario credibility)

6. Because the model runs on a personal computer within minutes, new geological, petrophysical, geophysical or drilling data can be incorporated as the well is drilled. It then provides cost effective, actionable information to the drilling managers such that they can proactively control their environment and make competent decisions about the future of their drilling program.

The RT adjustment of a prediction model requires re-formulation of the Specific Inversion Problem (SIP) introduced before in (Madatov et al., 2000; Madatov, Sereda, 2000a) for this special case. Below we discuss the approach to the RT prognosis updating and the necessary modifications required for the SIP definition. We also schematically describe the algorithm of the SIP solution in a real time mode. We address readers to our previous article (Madatov, Sereda, 2005b) for getting more details about the background theory and algorithm developed for the Specific Inverse Problem.

3.1. Confidence of the pre-drill prediction as a function of current drilling depth and an approach to its real time improvement

Above, in section 2 we have investigated the origins and the sources of the errors in a prediction model and related to it pre-drill overpressure prediction scenarios. It has been shown, that the instability cases in numerical forward modelling and model parameter interpolation can be checked and fixed a priori.

The main causes of potentially non-removable prediction errors are associated with the inherent non-uniqueness in the e-equivalent SIP solution and also with the tolerance uncertainties of the framework model convertible to the finite depth errors at the point of prediction. Still the maximal level of the pre-drill prediction scenario error is surely limited for all stable cases by calibration accuracy, which is in turn constrained by data stop criteria. The non-stable cases are detectable a priori during stability tests of the definitive prediction scenario to the tolerant error in a framework model (see Fig. 3).

Nevertheless, in the general case when the framework model at the target position is not certain enough the inherent non-uniqueness of the SIP solution and apparent modelling errors could jointly create an envelop of uncertainties in prediction scenario, which exceeds accuracy of the pre-drill calibration.

Now we have to consider both groups: inherent model parameter uncertainties and tolerant for the framework model depth errors together and more formally. This is important in order to estimate their direct integral effect of current prediction uncertainty as a function of current drill bit position. Note that the prediction uncertainty function is supposed to be inversely proportional to a confidence function of the pre-drill prediction. Our consideration will be made on the level general enough for understanding of the key idea of the RT data inversion method and sophisticated enough for getting approach to the corresponding SIP modification.

3.1.1. The uncertainty of a current prognosis

Let the prediction model {x} have been formed based on a calibrated framework HD model and the accuracy of the corresponding SIP solution for the available multi-well data set has been conditioned by the area common data stop criterion e. Let also the definitive pre-drill overpressure scenario have been estimated at the target well position of the prediction model in the form of the vector pA.

There is some inherent uncertainty in the pre-drill overpressure scenarios even in case when the model parameters' interpolation is absolutely precise and the forward pressure modelling is absolutely stable. This is due to the inherent local non-uniqueness of the prediction model Sx and due to the unavoidable uncertainty 8z in the subsurface geometry image of the framework model at the target position. The origins of these inherent

uncertainties in the prediction model and stability of the relevant prediction scenario at the target position were considered in details in section 3. Since all unstable cases can be controlled and prevented we will discuss here only a stable case, where we will recognise it as an inherent component of prediction error: Z ^ e.

The following comment must be premised at this context.

Despite different origins of the model vector uncertainty Sx and the framework model uncertainty Sz there is no needs to keep them separately in consideration of a stable case in pre-drill overpressure scenarios. Indeed, any top depth-thickness error attached to a certain formation unit can easily be converted into the relevant sedimentation rate error and then be compensated by the relevant addend(s) into sensitive component(s) of the formation control parameter(s) (Madatov, Sereda, 2000a). For example an apparent increase in sedimentation rate leads to the corresponding acceleration in pore water release during the rock matrix compaction. In order to compensate this effect on basin scale hydrodynamics of the given formation it is enough to increase its HD conductivity by tuning specific surface constant. The result in terms of a final pressure drop/build up won't change.

Due to the local continuity of the specific inverse operator inside of the Z-locality of pA it is easy to link it with some S-locality of the model vector prototype xA within the model parameter space (see Definitions 4-5 in Appendix 1). The uncertainty interval SxA, j=1,2,...M for every j-th component of the model parameter vector xA including the framework model uncertainty can be represented in the form:

SxAj = { xj : (xA- - SxA) < x< (xAj + SxA), VSxA- e {SxA-} }, j = 1,2,.. .M;

i

M = Z N,

i= 1

(7)

where xAj is the j-th component of most likely estimation of the model parameter vector xA achieved at the target position; I is the number of the formation units containing uncertain model parameters which require adjustments during subsequent drilling; Ni is the number of sensitive control parameters for the i-th formation unit.

Clearly, that the contribution of more uncertain component of the model parameter vector xA to the inherent overpressure prediction error Z is higher than it is for more certain component at the other equal conditions. This effect can be prior estimated approximately based on the local sensitivity analysis (Madatov, Sereda, 2003).

Note, that the magnitude of the j-th component of the sensitivity vector Sj computed according to formula (6) represents the total misfit between reference and current synthetic data curves averaged along the whole depth interval according to the vector norm definition (see Definitions 1 in Appendix 1). In the situation with the RT prognosis updating it is important to involve the Differential Sensitivity Analysis (DSA) instead of the standard or integral sensitivity analysis. This kind of sensitivity analysis allows quantitative estimation of the influence of the given control model parameters xj attributed to the given i-th formation unit on the synthetic pressure response as a function of the depth distance between this formation Zi and the current position Z of the DSA window. The changes in the sensitivity vector magnitude for certain j-th component and the fixed DSA window length can be represented as a continuous function: Sj ~ Dj (| Z, - Z |) of the current distance to the target formation | Z - Z | -Distance Influence Function (DIF). We address readers to Appendix 4 for getting the relevant definitions, formulas and illustrations. Here, we just stress that the noticeable DIF magnitudes correspond to the most sensitive x component(s) (key component(s) of the given HD sub-class). They are generally detectable for all sub-classes of the HD models faraway from the relevant formation unit along the DSA pattern. This forms prerequisites for usage of the pre-drill prediction model at the real time prognosis updating. Yet, the DIF gives a convenient and quite general approximation of the rigorous sensitivity vector calculation in the integrated form. This approximation appears to be helpful in particular for the general description of the current prediction credibility changes in a real time situation, when part of the model and real data uncertainties are converging to zero during drilling.

Since the DIF gives pressure response normalised to the standard model parameter variation (see Appendix 4) it can be used for estimation of the current integral prediction uncertainty (error Z) as a function of current drill bit position.

Let the target formation unit be marked with index "i" and have a top depth level relative to the reference surface equal to Zi. The contribution of the uncertainty in the particular xj model vector component related to the i-th formation (i j-th component of Sx) into the prediction uncertainty at the depth level Z can be approximately expressed via the relevant DIF component as the following:

Ap (xj, Z) * [W(xj) D(x; | Z - Zt |)], (8)

where D(xj | Z - Z{ |) is the i j-th component of the Distance Influence Function (see Appendix 4 ); W(xj) is the DSA window function proportional to its length.

Let us the drilling process in a planed target well have been started from the relative zero depth level and the current drill bit position (current RT moment at constant drilling speed) be equal to Zc>0. Formally it allows to substitute the pre-drill definitive overpressure scenario above this level with the corresponding real

pressure data (see Fig. 4). Thus, the prediction uncertainty does not formally exist above the current level ZC.

Below ZC the prediction uncertainty formally appears again2.

The dynamics of the uncertainty function AP) (x,, Z) as ZC increases can be expressed as follows:

0, 0 < Z < ZC

AP(x, Z) = ^ Mc Ni (9)

E E APi, (x,, Z), ZC < Z < Z*,

i=ij=\

where Z* and ZC are the planed maximal depth and current drill bit depth at the RT drilling moment respectively (Z* > ZC); Ni is number of sensitive control parameters for the i-th formation unit; Mc and M are the total and remained for penetration number of the formations respectively. Clearly, that Mc < M since at least one of the top formation units has already been penetrated as soon as ZC > 0.

Now let us express AP) (x, Z) explicitly in the form of a normalised error by using formulas (8-9) and approximation proposed for for D(x,; |Z - Z,|) (see formula A4.2 in Appendix 4 ). By this one can get the following expression for the current relative prognosis error estimation:

0, 0 < Z < ZC

E(x Z) = J mc m (10)

E E W(x,) S, e-[(Z - Zi)'Bi]/ pA(x,Z), ZC < Z < Z* ¿=1 j=i

where pA(x,Z) is the Z-grid component of vector pA representing the definitive prediction value at Z level.

Note that the total number of additive elements contributing in sum (i0) is decreasing as the current drill depth Z increases. If the sensitivity vector components are distributed along the formation uniformly then we can expect a monotonous decrease of the prognosis error at the given depth level Z as the drill bit will approach to the bottom level Z*. I.e. formally [E(x, Z)] ^ 0.

Z*

When the reference coordinate system is fixed not at the given depth level Z but at current drill bit position ZC this situation can be reviewed and represented in the form of the speculative cone of prediction uncertainty (Fig. 4). Here the gradual enlarging of the uncertainty cone with the depth is due to amount of the key overpressure factors controlling the present day overpressure regime grows up with the depth since the number of additives in sum (10) increases. Roughly this bunch of scenarios form starts the cone of prediction uncertainty (see Fig. 2e,f and Fig. 4a). As the drilled depth interval increases the amount of fixed framework model features and the number of certain control model parameters in the prediction model also increases. Thus, the cone of uncertainty will be gradually shrunken (compare Fig. 4c and 4d).

A sharp change in overpressure scenario can normally be caused by influence of some single dominant overpressure factors. The relevant control model parameter of the prediction model should have maximal sensitivity and longest DIF response (see Fig. A4.1 in Appendix 4 ). Hence it can be detected in advance at shallower depth as the drill bit will approach to the relevant formation.

3.1.2. A real time confidence function

As it follows from (10) the relative RT prognosis error in the monotonous case can be approximated by the Gauss-like function:

EX(Z) « (1/ AX) e(Z- Z*)2 ' Bx , ZC < Z < Z* (11)

where AX and BX are some adjustable parameters depending on the type of the calibrated HD model (lower index "X"); Z* and ZC - the planed maximal depth and current drill bit positions of the drilling well respectively (Z* > ZC).

Note, that the discontinuous case in the RT prediction error dynamics can be associated with a gross error in the framework model at the target well position such as missing of some formation(s), or missing of a significant unconformity(ies) within the predicted section interval, or misleading in litho-stratigraphy interpretation for existing formation units, or wrong tectonic compartmentalisation and so on.

Here we assume the errors in the framework model to be continuous and tolerance.

2 In practice the pressure measurements while drilling (MWD, RFT, LOT etc.) are only possible inside good permeable depth intervals associated with aquifer formations. Generally, post drill pore pressure estimations are based on multi-source data analysis, which is not trivial and obviously implies the interpretation factor. Hence it is not free from errors (see, for example, Mouchet, Mitchell, 1989; Doyle, 1998). The same is true for the framework model uncertainties, which could be only reduced during drilling, but not fully removed. Still, the problem of the post drill data analysis and accuracy of the posterior earth model reconstruction is far beyond the scope of the discussing problem. Here we just have to assume that the drilling process provides a geologist with absolute truth about prior predicted data.

1 1.25 1.5 1.75; 2 2.25;

"Pressure;[g/cc]

Fig. 4. Scheme of the while drilling prognosis updating in real time

a - the pre-drill overpressure scenarios with the definitive curve indicated (solid green line); b - RT drilling case 1; c,d - RT drilling case 2 before (c) and after (d) correction; e - the final fit between the real pressure data (black markers) and current overpressure scenarios. The dotted red line on (b-d) indicates the weight confidence function CX(Z) (see formula 12)

:

nU

i

Pressure, [Uff

b

1 75 2 2.25 Pressure, |||¡¡§

d

1 2: é;.2S HHpure

The normalised confidence function to the current prognosis can be therefore expressed in the form:

Cx(Z) = { A «/ If ZC < f . (12)

[ Ax/ Ex(Z), Zc < z f z ,

where all notifications have the same meaning like in formula (11).

Now the range of variation of CX(Z) belongs to the [0,1] interval. It is constructed in the form to be inversely depended on the current prognosis error (11) and to be maximal (CX(Z) = 1) at the level Z f Zc, where the real pressure data are already made available.

a

e

3.2. Modification of the Specific Inverse Problem for the while drilling purpose

The Specific Inverse Problem (SIP) has been introduced in our previous papers (Madatov et al., 1997; Madatov, Sereda, 2000a) as a non-linear data fitting problem (Lawson, Hanson, 1974) required at calibration of the framework model against single and multiple well data set(s). The approach and strategy of the SIP solution based on local optimisation principles recently were extended by global regularisation criteria and account for the HD model constrains (Madatov, Sereda, 2005b).

The RT prognosis updating can be treated as a local re-calibration problem in regard to the near well locality of the prediction model established before drilling. In agreement with the prediction strategy introduced above some certain HD model sub-class has to be predefined before calibration of the relevant framework model. Thus, its local re-calibration can be considered as a refinement of the conditionally correct local SIP solution against any new pressure related data, which can be made available during drilling.

Indeed, there are two possibilities at the RT data versus prognosis verification: real data vs. predicted data misfit is tolerant or is not tolerant to the predefined current prediction accuracy e. These scenarios can be expressed formally as the following:

p] ¿ e : current prognosis is valid, (13)

RvtpA, p] > e : current prognosis needs to be updated.

Here ] is normalised mismatch function given by

RN[pA, p*] = II M[zA;xA] - p* ||/||M[zA;xA]||, (14)

where pA = M[zA;xA] and p* are predicted and actual pressure related data computed and observed above current drill bit position ZC respectively.

The inverse problem in connection to the second option in scenario (13) can formally be posed as the following:

x: {Rx[M(x),p*] < e, V x e Xj}, (15)

where Xj is the compact sub-space corresponding to a certain HD model sub-class predefined for production of the given prediction model within the model parameter space X.

This definition is quite similar to the general SIP definition given in (Madatov, Sereda, 2005b) for calibration of the whole framework model, but now it has two specifically RT features, which could be treated as a particular case. They are the following:

1. The SIP is now posed not for the whole model parameter space X, but for its local sub-space XJ containing the current prediction model xA from series of available {XJ} c X.

2. The real data p* currently available for fitting (15) represents not whole, but only upper part of the section where Z < ZC.

The first feature simplifies the SIP definition since it is reduced now to the locally unique e-equivalent solution. This is because of the fact that XJ normally represents the simply-connected compact vector sub-space (see Definition 3, Appendix 1). The only wrong interpretation of calibration results or gross errors in the framework model at the target location can potentially force a predictor to change the choice of the predefined pre-drill prediction model XJ. It is important however to keep possible errors in the framework model distinguished from the wrong interpretation of calibration results. Indeed, the wrong formation top/thickness prediction at the target well location could be promptly corrected while drilling in the prediction model without re-calibration of the relevant framework model. If the RT interpretation of drilling result reveals rather big mismatch between the prediction scenario and the real data jointly with the unexpected geology section, then validity of the prediction model becomes questionable. In the worst case when the target locality of the prediction model belongs to a totally different tectonic or/and deposystem block, the corresponding pre-drill prediction scenario is inapplicable for further RT updating. Still, even in this situation the RT correction can be made quickly by substitution of the current local sub-space XJ containing the current prediction model xjA by more suitable prediction model xlA and its host sub-space Xl of available from series of predefined prediction models (see sub-section 2).

The second restriction can be extended to the more general SIP definition given in (Madatov, Sereda, 2005b) by using introduced above normalised confidence function (12) as a continuous weight function of the current depth Z. Indeed, the vector norm used in the mismatch function (14) can be extended to the weighted form. The weight coefficients can represent, for example, a confidence to the real data accuracy, which is normally not equal along well path. The higher confidence weight the more contribute the relevant component of data vector to the total misfit and vice versa.

Let the weight coefficient at the misfit function be assigned as a depth function Cx(Z), which has firstorder discontinuity3 in agreement with prediction confidence formula (12). Then the goal minimisation functional in (15) is getting the following RT form:

Z* Z*

R(x,Z) =!{Cx(Z) ■ [M(x) - p*(Z)]2} / X[Cx(Z) ■ p*(Z)]2 , (14*)

Zi=0 Zi=0

where the weight function Cx(Z) is getting by formula (15) and the forward modelling operator M(x) represents synthetic pressure data.

Evidently, all real data available at the current drill bit position above ZC depth level contribute to the total misfit function (16) with the equal weight (Cx(Z) =1 at 0 < Z < ZC; see also footnote 2). Below this level the real data are not available yet. Still, the most likely prediction vector components could formally be added as missed data components with a lower confidence. The contribution of predicted real data into the total misfit (14*) must quite sharply and non-linearly decrease with the depth (see Fig.4b-d). This provides both consistency of the RT modification with the general case and account for gradual increase in prediction uncertainty below drill bit (cone of uncertainty). It is important to pick out that goal function (14*) in formula (15) remains to be continuous, which insure numerical stability of optimisation process.

The strategy of the RT data inversion algorithm can be summarised as the following.

The current definitive prognosis ahead of the drill bit is getting by running of the forward overpressure modelling based on the current SIP solution. Since xA is defined for the whole depth interval of the drilling well, the relevant synthetic data pA = M(xA) will also be available for whole interval. At this:

3 Discontinuity in the first order derivative

The upper part of pA(z) at 0 < Z < ZC belongs to the e-locality of available data p (z).

The lower part of pA(z) at ZC < z < ZB forms the definitive prognosis below drill bit.

The uncertainty of prognosis ahead the drill bit is represented in Fig. 4a-d in the form of confidence intervals (light yellow areas) attached to the definitive scenario.

3.3. A robust solution of the real time prognosis updating problem

As it follows from the above definitions, the solution of SIP in the real time mode can be treated as a local optimisation problem (see also Definition 10-11 in Appendix 1) determined for the local sub-domain Xj predefined in the EBM model parameter space X. The main formal difference between the single well data inverse problem considered earlier (Madatov, Sereda, 2000a) and this modification is reduced to updating of the misfit function (14) to the form (14*). Thus, the core of the numerical algorithm remains the same, but the logic of internal processes' commutation is updating as the content and amount of the real data necessary for inversion assumed to be constantly enlarging.

It is important to pick out that the input information and data available in practice during drilling are represented in the form of geology section interpretation and petrophysical measurements (wire logging; core analysis; MWD, RFT, LOT etc.). The first group of this input (input information) allows recognition and adjustment (if necessary) of litho-stratigraphy prognosis (vector zA) at the target location. The second group of the input (input data) usually contains porosity and pressure related data applicable for verification of the current ahead of the bit prediction (vector xA) and its local re-calibration (if necessary). We address readers to the following references (Mouchet, Mitchell, 1989; Doyle, 1998) in order to get more detailed description of the second group data input.

The content and consequence of necessary data-model processing operations associated with while drilling RT control and updating of the definitive overpressure prediction can be expressed roughly in the form of following four steps computer algorithm:

STEP 1. Get new data applicable for adjustment and/or re-calibration of the current prediction model.

STEP 2. Check the current misfit between whole collected for the moment real data and their synthetic prototype computed basing on the current prediction model.

Drilling case 1: The current misfit is acceptable (the current prognosis is valid in (12)) ^ Stand by till new data portion comes from the well site. When it is happened - go to step 1. Drilling case 2: The current misfit exceeds prediction error threshold (the current prognosis is not valid in formula (12)) ^ Go to step 3.

STEP 3. Recalibrate the current prediction model in order to minimise the current misfit down to the acceptable level in agreement with the definition of the SIP solution for while drilling mode (formulas 13-15).

STEP 4. Use the updated synthetic overpressure curve below the current drill bit position as an updated version of the definitive prediction.

These steps are illustrated by the speculative example in Fig. 4. The RT prognosis updating stage is activating as soon as the very first during drilling data can be made available. The geology data allow verification of the current formation model for the target well position (vector z). Note that the seismic while drilling data potentially allow to perform this verification in real time ahead of the drill bit. The well data (measurements, logs, indicators etc.) convertible to porosity and pore pressure scales allow reconstruction of the definitive real data set and checking the current model vs. data misfit. There are potentially two possible cases in connection to these checks mentioned above as drilling case 1 and drilling case 2 represented in Fig. 4a-b and Fig. 4c-d respectively.

Let us reconsider them more formally.

Drilling case 1: The current misfit is acceptable.

Neither model consistency misfit nor real synthetic data misfit exceeds the predefined prediction quality stop criteria. In this case both check criteria are satisfied and current synthetic curves can be used as a prediction scenarios ahead of the bit (Fig. 4b-c). Thus, the ahead of the bit prognosis will be checked but not updated.

Drilling case 2: The current misfit exceeds the prediction error threshold.

At least one of the check criteria is not satisfied. In this case the current prediction model must be corrected and the relevant definitive overpressure scenario must be recomputed. Thus, the ahead of the bit prognosis will be updated basing on the fixed misfit (Fig. 4c-d). There are three sub-cases possible at this scenario:

Sub-case 2.1 (adjustment of the prediction model geometry): The initial formation model setting was changed in terms of lithology or/and stratigraphy or/and top formation depth based on the RT interpretation of the drilled section. The misfit between overpressure scenarios computed prior (definitive) and posterior (current)

these changes exceeds the predefined prediction quality stop criteria. If so, then the current local framework model z must be adjusted in agreement with the new formation model corrections.

Sub-case 2.2 (local re-calibration of the prediction model parameters): The misfit between real data and current prediction scenario exceeds the predefined prediction quality stop criteria and the initial formation model setting is in the tolerant agreement with the current drilling results. If so, then the current prediction model must be re-calibrated by using the current misfit as an input to the data inversion engine.

Sub-case 2.3 (adjustment of the prediction model geometry and local re-calibration of the prediction model parameters): The misfit between real data and current prediction scenario exceeds the predefined prediction quality stop criteria and the initial formation model setting was changed in terms of lithology or/and stratigraphy or/and top formation depth based on the RT interpretation of the drilled section. If so, then the current framework model must first be corrected according to sub-case 2.1 and then the current prediction model must be re-calibrated according to sub-case 2.2.

3.4. The speculative and real data tests of the method

We have used the test area with 12 calibration wells (see Fig. 5) in order to test the converging speed of the described above computer routine and also to estimate achievable quality of the RT prediction. Naturally, the real time process was artificially simulated for every real well by portion "opening" of the available information attached to the well. At this every new piece of the data including description of the over-drilled section, estimation of porosity and pore pressure were getting available in the format of the required input files. This input was continuously compared with the current definitive pre-drill prediction. Depending on the current misfit level both of the RT drilling cases described above were simulated.

The tests were made basing on the same North Sea calibration area and the real data multi-well data set, which has been used in our previous paper (Madatov, Sereda, 2005b). We address the readers to the following papers (Chiarelli, Duffaud, 1980; Doyle et al, 2003; Cubala et al, 2003) in order to get more details about the pore pressure data interpretation and geology setting within the region.

The summary of the well data used for the tests is available in table 1, where the names of the well correspond to the map in Fig. 5 and the name of the formations indicated in the parenthesis correspond to the formation column shown in Fig. 2a-b.

A blind test strategy has been implemented at every of calibration wells considered to be a target well during examination (Standifrd, Matthews, 2005). Namely, the well chosen as a target well was temporary removed from the framework model and the relevant prediction model was re-calibrated basing on the second sub-class of the HD model (see Fig. 2b,d,f and Fig. 3b) in order to ensure a stable pre-drill prediction. Then the SIP solution for the remained multi-well data set was used in order to get the relevant prediction model set. The model parameters were getting at the target well location based on the interpolation between the remained calibrated wells in agreement with the prediction strategy (see also Fig. 1b-c). We used the worst possible combination of the model parameters (maximal or minimal prediction scenario in Fig. 1d) belonged to the relevant compact sub-set of the possible local SIP solutions in order to compute definitive prediction scenario at the target position.

Table 1. The data input available for simulation of prediction model re-calibration in while drilling mode

Well Name Data available for definitive curve interpretation

Overpressure scale Porosity scale

Measurements Indirect data Measurements Related well logs

023 RFT (F2, F5, F8, F11) EMW, Kicks, Gas log (F2, F5, F8, F11) BCS

035 MWD CDE, Gas log BCS, FCD

101 RFT (F2, F5); MWD EMW (F2, F5) BCS

105 RFT (F2, F5, F8, F11) EMW, Gas log (F2, F5, F8, F11) BCS, CN

807 RFT (F2, F5); MWD CDE, Gas log (F2, F5) BCS

804 MWD CDE, Gas log BCS

202 RFT (F5) EMW (F5) BCS, FCD

801 RFT (F5, F8) EMW, CDE (F5, F8) BCS, FCD

802 RFT (F5, F8) EMW, CDE (F5, F8) BCS, FCD

803 RFT (F5, F8) EMW, CDE (F5, F8) BCS

401 RFT (F2, F5); MWD Gas log BCS

301 RFT (F5); MWD EMW, Kicks (F5) BCS

Abbreviation used in the table: RFT - Repeated Formation Test; MWD - Measurements While Drilling; EMW - Equivalent Mud Weight; CDE - Compensated Drilling Exponent; BCS - Borehole Compensated Sonic log; CN - Compensated Neutron log; FCD - Formation Compensated Density log.

^^ Error in definitive pre-drill overpressure prediction; ^^ Error in real time prognosis correction; — Depth level of RT correction.

Fig. 5. Results of the speculative tests of the computer routine for the real time correction of the overpressure prediction in while drilling mode

The RT correction plots represent current prediction error (the absolute value real data - prediction misfit) versus drilling depth. The final RKB bottom marked with the grey bar. The correction occurrences marked with pink markers. The position of represented wells on the area map is highlighted with the stars. The rest of the tested wells are not presented on plots. (See full test summary on tables 1-3)

Table 2. The multi-well test summary

N Well name Number of the required RT corrections The depth of the last available data [m] The depth of the first correction [m] The depth of the last correction [m]

1 105 3 3855 1030 3210

2 202 2 3982 2250 3180

3 035 3 4508 3030 3700

4 023 4 4664 2950 4020

5 101 4 4783 2970 4190

6 301 2 4185 2020 3000

7 401 2 4577 2360 3720

8 804 1 3455 1900 1900

9 803 3 3923 1280 3900

10 802 4 4162 1790 3650

11 801 2 3942 1170 3850

12 807 3 5215 3070 4810

Three drilling sub-cases described above have been analysed:

1) error in the prediction model geometry and RT adjustment of the formation top(s);

2) error in the prediction model parameters and RT re-calibration;

3) combination of 1 and 2.

The results of the tests are summarised in table 2, Fig. 5 and in the below summary. • The pre-drill prediction can be effectively updated during drilling on condition that the interpretation of calibration results remains valid at the target position and the only local re-calibrations and/or adjustment of the

framework model are required. Otherwise, the updating (adaptation) of the definitive prediction scenario ahead of the drill bit could be significantly delayed. Thus, the more consistent choice of the prediction model sub-class to the real deposystem setting the more adaptive in real time the relevant prediction scenario.

• The adaptability of the pre-drill prediction can naturally be estimated in terms of the total number of the required prognosis corrections. This characteristic and averaged time used for the single correction can be recognised as most important indicators of the prediction quality. At this the total number of RT corrections is more controlled by the SIP solution HD sub-class predefined for the prediction model at other equal conditions.

• The tolerance errors in the local setting of the framework model geometry or/and lithology (drilling case 1) can distort prediction longer and can require more corrections than locally wrong parameterisation of the locally correct framework model (drilling case 2).

• Practically the most likely is drilling case 3, that combines two previous ones. Sill, the converging speed of the RT correction routine is high enough to produce updated version of the definitive prediction ahead of the drill at any realistic drilling speed

• The tests performed for all available calibration well reveal full scale RT adaptability for any of prediction scenarios produced basing on the globally consistent prediction model. In particular, the following adaptability estimations have been obtained (see also Fig. 5).

Table 3. The prediction adaptability summary based on the tests

No 1 2 3 4

Number of corrections required for whole well path 4 3 2 1

Number of wells / total number of wells tested 3 / 12 4 / 12 4 / 12 1 / 12

The very first real data examination of the described above method and computer procedure was made in 1997 during drilling of a real exploration well in offshore of the North Sea. The successful results achieved there allowed saving one casing in the drilled well, that in turn significantly reduced the total drilling cost and time. These results were briefly described in (Madatov et al, 1998). More general real data test of this method was made recently basing on deepwater wells in the Gulf of Mexico. Here the quality of Look Ahead predictions based on the described approach was considered in comparison with the standard method. The blind tests were performed over the entire well trajectory and involved more than 30 drilled exploration wells. The results of these blind tests were that 62 % of the wells (22 of 32) reveal accuracy of the RT prediction within ±1 ppg4 over the entire well path. At this, 50 % (16 of 32 tested wells) reveal accuracy of the RT prediction ±0.5 ppg for the same conditions. This accuracy significantly exceeds achievable accuracy by standard seismically driven method. More detailed description of the results is available in (Standifird, Matthews, 2005).

3.5. Discussion

Formally, the real time solution of the SIP in while drilling mode provides quick converging to the optimal look ahead prediction. However, if the choice of the background prediction model x for the target well position from the list of available {x} after calibration of the framework model is not optimal (see prediction case 1 in subsection 2) then the speed and quality of the RT correction won't be optimal too. The reason for this is flexibility of the EBM during its tuning against real data and finite accuracy of the while drilling data available also non-regularly. The wrong choice of the HD model subclass during recalibration of the current prediction model will be in this case compensated by production of the specific solution on the target position xA, which will be less and less consistent with the calibration results for the offset wells {xA} as the drilling process will proceed. Thus, it is important to submit that the adequate choice of the calibration area at producing prediction model becomes a paramount and the only prerequisite for success of the whole prediction strategy at the final while drilling stage.

We suggest to keep the local errors in the framework model (drilling case 1) or in the prediction model (drilling case 2) distinguished from wrong interpretation of calibration results. Indeed, the wrong formation top/thickness prediction at the target well location as well as the worst choice among e-equivalent local solutions could promptly be corrected while drilling in the prediction model without full re-calibration of the relevant 3D framework model. If the interpretation of drilling result reveals rather big mismatch between the prediction scenario and the real data jointly with the unexpected geology section, then validity of an applied prediction model becomes questionable and full scale re-setting and re-calibration of the basis 3D framework model could be the only helpful issue. In the worst case scenario (prediction case 1 in sub-section 2) the target locality of the prediction model belongs to a totally different tectonic or/and deposystem block. It means that region chosen as a calibration area is inconsistent with the target drill location. Hence, the corresponding prediction model sub-type is inapplicable for producing adaptable in real time prediction scenario.

4 ppg - pound per gallon is unit of a pressure gradient. The conversion formula is the following: 1g/cm3 = 8.33 ppg

4. Conclusion

The considered approach to the pre-drill prediction and its real time correction during drilling allow to extract major benefits of the effective basin model concept (Madatov, Sereda, 2005a) and multi-well data inversion strategy (Madatov, Sereda, 2005b) at the final stage of generic prediction strategy - during drilling prognosis updating. These benefits were emphasised before in connection with easier and faster 3D HD modelling and practically workable data inversion routine respectively. However, the only real time implementation of the developed approaches to the forward and inverse problem solutions allows to realise the full-scale advantage from usage of fit-to-purpose background models and concepts.

The results in the context of discussed real time implementation can be briefly summarised as the following:

• The prediction strategy is formulated basing on implementation of locally unique calibration results to the common framework earth model. The set of prediction models allows to produce the relevant set of the pre-drill pore pressure predictions (overpressure scenarios) based on combine 3D HD forward problem solutions. At this, each particular prediction model represents a local e-equivalent solution achieved for the framework area based on multi-well data inversion for the specific model sub-class Xj. The given prediction model then can be recognised within the model parameter space X by its unique combination of the highly sensitive and low sensitive interval of control formation parameters.

• Each of the produced prediction models can be evaluated basing on the formal confidence and stability criteria and non-formal criteria of its consistency with the deposystem analysis. The best of these model serves as a background for calculation of the definitive pre-drill prediction. Still it is highly recommended to keep all available prediction models ready for substitution of the current one during drilling correction.

• The RT modification of the Specific Inverse Problem posed for definitive porosity and pressure data implies determination of data confidence weights as a function of the current drill bit position. This function has the 1-st order discontinuity tied to the deepening while drilling drill bit depth level. Namely: the weight function has a constant level above drill bit position and exponentially decreases below it (see Fig. 4b-d). Relevantly corrected goal mismatch function of the SIP minimisation problem allows implementation of the developed earlier (Madatov, 2000a) approaches and algorithms for getting conditionally correct and locally unique inverse problem solution.

• The strategy of the real time updating of current overpressure prediction implies consequent application of the SIP solver to the pressure related well data routinely available during drilling along available depth interval. Any disagreements between formation model geometry or lithology and the real geology data obtained during drilling shall also be routinely entered into the prediction model in the form of the corresponding corrections.

• The speculative and real data tests revealed the real time adaptability for any of prediction scenarios produced based on the globally consistent prediction model. At this, the converging speed of the relevant computer program is high enough to produce updated version of the definitive prediction ahead of the drill at any realistic drilling speed

Basing on the achieved results we can state that the time spent on earth model setting, upscaling and calibration is returning back during drilling, when it costs much higher.

List of references

Braverman E.M. Computer tests in pattern recognition. Automatics and Telemechanics, v.XXIII, p. 13-21, 1962 (in Russian).

Brekke K., Thompson L.G. Horizontal well productivity and risk assessment. Paper presented at the SPE

Annual Technical Conference and Exhibition held in Denver, Colorado, USA, 6-9 October, 1996. Brevik I. Rock model based inversion of saturation and pressure changes from time-lapse seismic data. 69th

Annual International SEG Meeting, p.1044-1047, 1999. Cheney W., Kincaid D. Numerical mathematics and computing Brooks. Cole Publishing Company. Monterey,

California, USA, p.562, 1985. Chiarelli A., Duffaud F. Pressure origin and distribution in Jurassic of Viking basin (United Kingdom -

Norway). The American Assoc. of Petroleum Geologists Bull., v.64, N 8, p.1245-1266, 1980. Cubala M., Bastow M., Thompson S., Scotchman I., Oygard K. Millennium Atlas: Petroleum geology of the

central and northern North Sea. 2003. Dennis J.E., Scnabel R.B. Numerical methods for unconstrained optimization and nonlinear equations.

Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 440 p., 1983. Doligez B., Bessis F., Burrus J., Ungerer P., Chenet P.Y. Integrated numerical simulation of the sedimentation heat transfer, hydrocarbon formation and fluid migration in a sedimentary basin: The Temispack model. In "Thermal Modelling in Sedimentary Basins" collection of papers from 1-st Exploration Research Conference. IFP, France, p.149-172, 1985.

Doyle E.F. Data and model integration results in accurate quantification of geopressure at the wellsite. Paper presented on conference "Pressure Regimes in Sedimentary Basins and Their Prediction", Sept. 2-4, DelLago, Texas, USA, 1998.

Doyle E.F., Berry J.R., McCormack N.J. Plan for surprises: Pore pressure challenges during of a deepwater exploration well in mid-winter in Norway. Paper 79848 presented on SPE/IADC Drilling Conference held in Amsterdam, The Netherlands, 19-21 February, 2003.

Duppenbecker S.J., Iliffe J.E., Osborne M.J. The role of multi-dimensional basin modelling in integrated pre-drill pressure prediction. Paper presented on 11th meeting of the Geologic Modelling Society, Galveston, 5 September, 2002.

Dutta N.C. Geopressure prediction using seismic data: Current status and the road ahead. Geophysics, v.67, N 6, p.2012-2041, 2002.

Eide A.L., Omre H., Ursin B. Prediction of reservoir variables based on seismic data and well observations. Jour. Am. Stat. Assoc, v.97, N 457, p.18-28, 2002.

Gandi A.F., Berg R.R. Primary migration by oil generation fracturing in low permeable source rocks. The American Assoc. of Petroleum Geologists Bull., v.81, N 3, p.424-443, 1997.

Hadamard I. Le probleme de Cauchy et les equations aus derivees partielles lineaires hyperboliques. Hermann edition, 1932.

Huffman A.R. The future of pore pressure prediction using geophysical methods. The Leading Edge, February, p.199-205, 2002.

Hunter R.L., Mann C.J. (eds.) Techniques for determining probabilities of geologic events and processes. Oxford University Press, New York, 394 p., 1992.

Ivanov V.K. To linear ill-posed problems. RAS USSR, v.196, N 145, p.270-272, 1970 (in Russian).

Jackson D.D. The use of a priori data to resolve non-uniqueness in linear inversion. Geophys. J.R. Astr. Soc., N 57, p.137-157, 1979.

Landro M. Discrimination between pressure and fluid saturation changes from time laps seismic data. Geophysics, v.66, p.836-844, 2001.

Lawson Ch.L., Hanson R.J. Solving least squares problems. Prentice-Hall Inc., Englewood Cliffs, New Jersey, p.263, 1974.

Lerche I. Basin analysis, quantitative methods. Academic Press, San Diego, California, v.1, 562 p., 1990.

Madatov A.G. The overpressure driven seismic velocity response. Review of standard models and methods for extraction in the context of basin modelling approach to overpressure prediction. Proceedings of the Murmansk State Techn. University, v.8, N 1, p.84-119, 2005.

Madatov A.G., Sereda A.-V.I. A multi-well data inversion purpose-built for calibration of an effective basin model at overpressure prediction. Proceedings of the Murmansk State Techn. Univ., v.8, N 1, p.44-83, 2005b.

Madatov A.G., Sereda A.-V.I. The decomposition of 3-D overpressure evolution model in basin scale and its application to the fault seal analysis. Proceedings of the Murmansk State Techn. Univ., v.4, N 1, p.79-96, 2001.

Madatov A.G., Sereda A.-V.I. The effective basin model concept and fast 3D overpressure modelling in basin time scale. Proceedings of the Murmansk State Techn. Univ., v.8, N 1, p.5-43, 2005a.

Madatov A.G., Sereda A.-V.I. The forward and inverse problems of the fluid dynamics in basin modelling applied to the pore pressure prediction within the sedimentary basins. P.2. The practical aspects. Proceedings of the Murmansk State Techn. Univ., v.3, N 2, p.351-366, 2000b.

Madatov A.G., Sereda A.-V.I. Upscaling of a basin model required for pore pressure prediction before and during drilling. Proceedings of the Murmansk State Techn. Univ., v.6, N 1, p.119-144, 2003.

Madatov A.G., Sereda A.-V.I. The forward and inverse problems of the fluid dynamics in basin modelling applied to the pore pressure prediction within the sedimentary basins. P.1. Theory aspects. Proceedings of the Murmansk State Techn. Univ., v.3, N 1, p.89-114, 2000a.

Madatov A.G., Sereda A.-V.I., Doyle E.F. Integration and inversion of well data into the basin evolution model: Way to the new generation of pressure prediction technologies. Paper presented on forum "Pressure Regimes in Sedimentary Basins and Their Prediction" Houston, Texas, USA, 2-4 September, 1998.

Madatov A.G., Sereda A.-V.I., Doyle E.F. Pore pressure prediction by using inversion before and during drilling. Paper presented on Workshop "New Methods and Technologies in Petroleum Geology, Drilling and Reservoir Engineering", 19-20 June 1997, Krakow, Poland, 1997.

Magara K. Compaction and fluid migration. Elsevier Scientific Publishing Company, p.319, 1978.

Menke W. Geophysical data analysis. Discrete inverse theory. Academic press, New York, 312 p., 1984.

Merriam D.F., Davis J.C. Geologic modeling and simulation, sedimentary systems. Kluwer Academic/Plenum Publishers, New York, 352 p., 2001.

Miall A.D. Principles of sedimentary basin analysis. Springer-Verlag, Berlin, p.616, 2000.

Mouchet J.P., Mitchell A. Abnormal pressures while drilling, Manuels techniques. Elf Aquitaine, Boussens, 286 p., 1989.

Nelson G.M., Franke R. Surface construction based upon triangulations surfaces in CACD. North-Holland Publishing Company, p.162-177, 1983.

Nordgard Bolas H.M., Hermanrud Ch., Gunn Teige M.G. Origin of overpressures in shales: Constrains from basin modelling. The American Assoc. of Petroleum Geologists Bull., v.88, N 2, p.193-211, 2004.

Osborne M., Swabrick R.E. Mechanisms for generation of overpressure in sedimentary basins: A reevaluation. The American Assoc. of Petroleum Geologists Bull., v.81, N 5, p.1023-1041, 1997.

Slingerland R. Predictability and chaos in quantitative dynamic stratigraphy. In: T.A. Cross (ed.): Quntitative dynamic stratigraphy. Prentice-Hall, Englewood Cliffs, New Jersey, 625 p., 1990.

Standifird W., Matthews M. Real time basin modelling: Improving geopressure and earth stress prediction. SPE 96464 Paper presented at Offshore Europe Conference of Society of Petroleum Engineers held in Aberdeen, Scottland, UK, 6-9 September, 2005.

Swabrick R.E., Huffman A.R., Bowers G.L. Pressure regimes in sedimentary basins and their prediction. History of AADE forum. The Leading Edge, April, 1999.

Tarantola A. Inverse problem theory: Methods for data fitting and model parameter estimation. Elsevier, Netherlands, p.386, 1987.

Tikhonov A.N. About stability of inverse problem. RAS USSR, v.39, N 5, p.195-198, 1943 (in Russian).

Tikhonov A.N., Arsenin V.Ja. The methods of solution of ill-posed problems. Moscow, Nauka, 285 p., 1979 (in Russian).

Tissot B.P., Welte D.H. Petroleum formation and occurrence. Springer-Verlag, New York, p.538, 1978.

Traub J.F., Woznjakovski H. Theory of optimal algorithms. Academic press, New York, 183 p., 1980.

Ungerer P., Burrus J., Doligez B., Chenet P.Y., Bessis F. Basin evaluation by two dimensional modelling of heat transfer, fluid flow, hydrocarbon generation and migration. The American Assoc. of Petroleum Geologists Bull., v.74, p.309-335, 1990.

Vasco D.W., Datta-Gupta A., Behrens R., Condon P., Rickett. Seismic imaging of reservoir flow properties: Time-laps amplitude changes. Geophysics, v.69, N 6, p.1425-1442, 2004.

Wang H.F., Anderson M.P. Introduction to groundwater modelling. Finite difference and finite element methods. Academic Press Inc., 237 p., 1982.

Williams K.E., Madatov A.G. Analysis of pore pressure compartments in extensional basins. GCS-SEPM 25th Annual Research Conf., Petroleum Systems of Divergent Continental Margins, in press, 2005.

Xiaorong L., Vasseur G. Geopressuring mechanism of organic matter cracking: Numerical modelling. The American Assoc. of Petroleum Geologists Bull., v.80, N 6, p.856-874, 1996.

APPENDIX 1. Basic terminology A1.1. The important math definitions

The below definitions are in agreement with common math terminology (Hadamard, 1932; Tikhonov, Arsenin, 1979; Traub, Wozniakowski, 1980; Dennis, Scnabel, 1983; Menke, 1984; Cheney, Kincaid, 1985; Tarantola, 1987).

Definition 1

The N-dimensional vector space X is defined as metrical when p-metric as for its elements x_ e X is introduced as the following:

||x|| = (x/ + x2p + x3p +...+ x/)1/p, (A1.1)

where p - the index of metric. In particular, the Euclidean vector space has metric index p = 2 (L2-metric) and L2

Vector Norm

||x|| = (x12 + x22 + x32 +...+ Xn2)1/2. (A1.2)

We will only consider the L2 metric which, for scalar product, gives:

kvx2 )a = CX-1 Xi, (A

where the covariance matrix C has a diagonal structure, meaning that there are no linear relationships between the arbitrary vectors in any space.

Definition 2

The set AX c X is defined as a Compacted Set in X if any sequence of its elements {x} contains subsequence converging to single element x0 e X.

The Diameter (size) of compacted set AX is defined as:

0(AX) = sup (Hx- - xj ||) (A1.4)

Vxj,XiC8X

Definition 3

The Compacted Set {x} in X space is defined as a simply-connected (or self-consistent) in X if it is closed, i.e. if any sequence of its elements {x} is converging to single element x0eSX. This definition can be formally expressed as the following:

0< e1< s2< s3 : 0<S1<S2<S3 ^..0 * {S 1x} c {S 2x} c {S 3x} c Xu, (A1.5)

where 0 denotes Diameter (size) of compact set {x} defined as:

0 (SX) = max (Hx - £ ||) (A1.4*)

Vxj,xjcSX

Definition 4

An s-locality {ep0} of the arbitrary vector p0 e P is defined as compact eP c P centered to p0 such that:

eP| p0 = {ep0}: max(ta- sj J< e. (A1.6)

vD.ceP V y

U-i

Definition 5

Let us there be two Euclidean spaces X, P and on some sub-space X0 c X the rule of imaging G: X^P is defined to map X0 into some sub-space P0 c P. The G[x] will be specified as an Mapping Operator (or just Operator) such that:

V x e X0 ^G[x] = p e P0. (A1.7)

Below we will assume that Vp e P0 3 x e X0.

Thus, will P0 c P call Image of X0 c X on the vector space P and X0 - Prototype of P0 on vector space X. Definition 6

The operator G[x] performs a Single-valued (Unique) mapping of the vector x e X0 if the unique image G[x] = p e P0 exists.

Definition 7

The operator G[x] performs a Continuous mapping of the compact set SX c X0 on P0 c P if reducing the size of SX (see (A1.5)) implies the relevant reducing of its image ePc P0. More strictly:

Sj< S2 = diam(SjX) < diam(SjX) ^ G[SjX ] c G[S2X ] ^ ejP c e2P. (A1.8)

Definition 8

The operator G[x] performs a Stable mapping of the vector x e X0 on P0 c P if small perturbations of x do not create big perturbation in its image p. More strictly:

G[x + Sx] = p + ep : ||SX || ^ 0 ^ ||ep || ^ 0. (A1.9)

The continuous on X0 c X operator G[x] is stable on the same sub-space.

Definition 9

The sub-set of all possible s-equivalent solutions {x}N of an inverse problem is representable as an aggregation of all pairwise non-crossing locally unique s-solutions given in the form:

{x}N = [{Sx}ju{Sx}2 u {Sx}3^u {Sx}n], (A1.10)

where {Sx}j * 0, j = 1,2,_,N; [{Sx},o{Sx}j] = 0, V i *j:i, j e (1,2,...N) is the diameter of compact sub-set (see Definition 3, formula A1.4*).

Then the condition of existence of global non-uniqueness in an inverse problem solution can be expressed as the following:

{3x c {x}N: xg{Sx}j Vj = 1,2,...N}. (A1.11)

Definition 10

The operator G[x] in this context is defined as Forward operator of bijective transformation (or just Forward operator) and operator G_1[pJ - as an Inverse operator of bijective transformation (or just Inverse operator). Both operators are Monotonous on the relevant sub-space of their monotony.

The wide class of linear (matrix) operators are having an inverse matrix and are monotonous wherever on the vector spaces of their definitions. The pair of bijective transformations can be written as the following:

{Gx = p e P - Forward transformation;

G_1p = xe X - Inverse transformation.

The typical example of bijective transformations and the relevant numerical operators are forward and inverse Fourier transform (Cheney, Kincaid, 1985). Generally, however, the unique and stable inverse operator G"1 [p] cannot be globally (wherever) defined on the Euclidean vector space P based on globally stable and unique forward operator G[x] defined on Euclidean vector space X. The relevant inverse problem is ill-posed in Hadamards' sense (Hadamard, 1932).

According to (Hadamard, 1932), the math problem Gx = p for p e P and x e X is well-posed (correct) if three following conditions could be accomplished:

{The solution can be find for the whole vector space P.

The solution is unique.

The solution is stable.

Definition 11

According to Tikhonov (1943) the inverse operator G"1 can be treated as a Regular Operator on subspace P0 c P if it is possible to build there such a stable G~l operator that the biggest diameter of the compact {Sx0} getting on the vector space X0 c X as a result of its application to compact {ep0} defined wherever on P0 is continuously decreasing with decreasing of e. The condition of regularity for solution of the relevant inverse problem defined on the sub-space P0 c P for the target sub-space X0c X can be then written as the following:

lim|| G"1e[pe] - x ||= 0,V pee{ep0}c P0; V xs e{Sx0}c X0. (A1.12)

When these conditions can be accomplished, the ill-posed inverse problem becomes Conditionally Correct.

The inverse problem can be treated as a conditionally correct according to Tikhonov if three following conditions can be satisfied:

{The solution of IP exists and belongs to a priori known sub-space X0 c X of possible solutions.

The forward operator G[x] is bijective for Vx e X0 c X; and G [x] = pA e P0 c P Vx e X0 c X.

The inverse operator is continuous on P0 c P.

A1.2. The important topic definitions

The Framework Model implies litho-stratigraphy unification and 3D subsurface architecture image established within the calibration area on the formation level of sequence hierarchy based on results of the source earth model upscaling. Namely:

lithology dependent porosity-permeability relationships defined for every formation units; unique stratigraphy column with all formation intervals' presented; present day top formation geometry image;

set of eustatic paleo sea level, paleobathimetry and temperature constants based on a priori given basin analysis results (see Madatov, Sereda, 2005a for more details).

Control (tuneable) HD model parameters are represented by the set of sensitive to overpressure parameters of the HD basin model adjustable during calibration inside their variation range by using a data inversion routine. Namely:

matrix compaction constants; HD conduction constant; HC generation constant;

pressure communication constants (for aquifer formation only).

The relevant model parameter space X accommodates full variability range for all of these control parameters (see Madatov, Sereda, 2000a,b for more details).

A sub-class of the fit-to-puaposes HD basin models is defined as a type of the physically possible HD basin model representable within the model parameter space X by compact set (cluster Xj c X) of the overpressure sensitive control parameters. Every specific sub-class XJ typically has its unique combination of the key overpressure controlling factors given by sensitivity vector components.

A Framework Model arranged with the sets of tuneable HD basin model parameters calibrated against of a multi-well data set and ready for computing of the relevant synthetic overpressure curve wherever inside the calibration area represents a Prediction Model. The Prediction Model must belong to one of sub-class of relevant for purposes HD basin models priorily accepted as a physically possible origin of the observed phenomena.

APredictionjicenaria

The single bunch of synthetic overpressure curves computed at the target well position and arranged with error bars forms a single prediction scenario (see Fig. 1d and Fig. 2e-f). This scenario is based on a full set

of locally unique e-equivalent inverse problem solutions associated with a certain Prediction Model.

The most lOvely peediction scenario uniquely defined basing on stability and consistency analysis performed over all a priori accepted Prediction Models (see Fig. 1d and Fig. 2e-f).

APPENDIX 2. The numerical instability at the model parameter interpolation

Let the framework model in agreement with the effective basin model architecture (Madatov, Sereda, 2005a) consist in set of surfaces reproducing geometry and model parameter variations attached to every formation unit and numerically represented in the form of stack of triangulation grids like it is shown on slide. Proposed for use interpolation function (Nelson, Franke, 1983) in Delone definition is given by:

r M(x,y) = C (x,y) + D(x,y)

2 2 C (x,y) = bi [Fi (3 - 2bi) + Ski bk - Sji bj]+bj Fj (3 - 2b)

^ + Sy bi - Skj bk] + bk2 [Fk (3 - 2bk) + Sjk bj - Sik bi] (A2.1)

D(x,y) = w {bi [Fi 6 (bkay + ^a^ + Ski F - ^ FJ + bj [Fj 6 (bi^k + bkaji) + Sy Fjk - Sk j + bk [Fk 6 (bjaki + biay) + Sjk FH - Sik F^]}.

where bp bj, bk - barycentric coordinates of the current point inside the arbitrary mesh forming triangle T(IJK) displayed in Fig. A2.1; w = b1bjbk / (bibj + btbk + bjbk);

"Sy = xjk Fjx+ yjk Fjy S = xki Fix + yki Fiy f Fj = 3bkaij + bj - bk

Sjk = ^k Fkx+ yjk Fky 4 Sji = xy FiX + yij Fiy J Fik = j + bk - bj

Sik = xki Fkx + yki Fky I Sy = xij Fjx + yij Fjy

Fjk = 3bajk + bk - bi

, = xj - xk; y,k=y - yk; xki = xk - xi; yki = yk - xij = xi - xj; yij = yi - y,;

aij = (ei + e, - ek) / 2e, aik = (ei + ek - e,) / 2ek a,k = (e, + ek - ei) / 2ek

aji = (e, + ei - ek) / 2e aki = (ek + ei - e,) / 2e akj = (ek + e, - ei) / 2e

F,i = 3bka,i + bi - bk Fki = 3b,aki + bi - b, Fk, = 3biak, + b, - bi

2 t 2 2 t 2 2,2 ei = x,k + y,k; e, = xki + yki; ek = x, + yi, •

This function provides a minimal derivative norm and optimal edge smoothness conditions for interpolating surface.

Still the optimal smoothness and minimal norm conditions can be violated inside an extreme slime triangle, if such one exists. The relevant non-stable constrained triangulation cases could appear when a target mesh knot is very close to one of existing knot or/and edge of the triangulation area

The mesh forming triangle (see Fig. A2.1) was recognised as an extremely slime one in case if its biggest side does not exceed the smallest one in seven times. I.e. when:

| JK | < 7 | IK |.

In this case an updated version of Nelson & Franke (1980) algorithm is suggested for implementation. This version implies the following substitutions in formula (A2.1):

K

Fig. A2.1. A mesh forming triangle

âji = 0.5; aki = 0.5; âj = (a^aji)/ âj; âik = (aki aik)/ âki; âjk = 1- âik; âkj = 1— âij,

This version ensures stabilisation of interpolation error like it shown in Fig. A2.2.

(A2.2)

a b

Fig. A2.2. Stabilization of the interpolation procedure on a triangulation mesh a - the raw interpolation of an irregular discrete field in continuous surface form; b - the same irregular discrete field interpolated by using new algorithm in continuous surface form

APPENDIX 3. The numerical instability at the final difference approximation of the forward modelling operator

This topic was in details considered in our recent article (Madatov, Sereda, 2005a). Here we just remind most important statements.

The Pieceman & Rachford scheme of the ADI method implemented in the forward modelling operator M [z;x] provides an absolute modelling stability in stable mapping sense (see Definitions 7,8 in Appendix 1) for every current solution achieved on a fixed time-spatial grid (Wang, Anderson, 1982). The numerical stability of the forward modelling operator M[z;x] in contrast is property of its final difference approximation, which ensures repeatability of synthetic data computations when simple steps of the grid are changing from run to run. As soon as vs. grid parameters can be treated as a stationary process that never exceeds certain modelling error level "e, the relevant operator M[z;x] determined jointly with its grid property can be treated as numerically stable.

In (Madatov, Sereda, 2005a) we have defined a discreteness threshold as maximal time-spatial sampling above which the numerical stability could fail. This inherent numerical modelling stability condition can be ensured wherever on the model parameter space X if the time-spatial samples do not exceed discreteness thresholds.

The results of the relevant numerical testing could be summarised as the following.

The combined 3D scheme implemented in the forward modelling operator M[z;x] provides reliable converging of current synthetic data to the reference one when the number of single charging-discharging cycles per formation burying with the maximum sedimentation rate 1000 m/M.Years is bigger or equal to 10. At that, the misfit value between synthetic data computed on such gross grid and the reference prototype is still below 1 % in average. Note, that this is much less than the generally accepted data stop criterion.

The pressure and porosity solutions achieved on less dense time grid are not surely stable. The behaviour of inherent modelling error as a function of the model vector position is not stationary. Depending on the model parameter range estimation of can randomly exceed 10 % level.

The single time cycle threshold for the sedimentation rate ranged between 50-500 m/M.Years is about 0.1 million years. Below this level the modelling error does not exceed 0.5 % and converging process becomes stationary. Above this level the averaged error gradually decreases.

APPENDIX 4. The Differential Sensitivity Analysis and Distance Influence Function

Let any uncertainty in the framework model at the given target position within the area be fixed or could be firmly compensated in a synthetic data by tuning of adjustable model parameter vector x. Than the sensitivity vector defined at any arbitrary position of the relevant EBM model sub-class e Xj c X can be determined as follows:

S(x) = ||MIx+dx] - M[x]||/||M[x] ||, Vx e XJ;

(A4.1)

where dx is infinitesimal local perturbation of vector x, such that ||dx||"||x|; M[x] = p e P - forward modelling operator; the L2 norm of synthetic data vector p e P is given in agreement with the dimension of data vector space P.

Let the sensitivity of the combine 3D HD model (Madatov, Sereda, 2005a) is investigating along single well path. Without changes in commonness of discussion we can assume for simplicity that the relevant well grid has only one dimension coinciding with the depth axis (vertical well case).

If the depth interval used for the vector norm computation in (A4.1) is fixed and is moving along the depth axis Z with the constant overlap step 8Z the sensitivity vector (A4.1) becomes depending on the current position of the norm interval - S(x,Z).

Let this modification of the sensitivity analysis be called the Differential Sensitivity Analysis (DSA) and the scanning along the Z-axis norm interval of misfit computation be called a DSA window. This kind of sensitivity analysis certainly allows closer investigation of the dynamics of sensitivity component redistribution for the specific target formation as the DSA window approaches to it (see Fig. A4.1). It also emphasises a mutual dependence between formation elements included into the common HD model. In particular it proves existence of noticeable impact of a remote formation unit on overpressure regime observed against certain depth level.

The computer tests of the Differential Sensitivity Analysis for the different sub-class of HD models with the big degree of confidence showed the following common features of a DSA signature:

1. The more sensitive j-th component of the vector x to the synthetic overpressure data p, the bigger depth interval of its noticeable influence in terms of the DSA magnitude.

2. The changes in the DSA magnitude for certain j-th component at the fixed DSA window length can be represented as a continuous function of current distance to i-th target formation - a Distance Influence Function (DIF) (see Fig. A4.1).

A Gauss-like function type is suggested to be used for approximation of the DIF in the following form:

E ■

E ■ ■ ■ ■ ■ ■ • 11 • * * • •

E • • • • * * * * * \

b

1.0 1.25 1.5 1.75

Sensitivity

2.0

Fig. A4.1. The differential sensitivity analysis of pore pressure model

a - synthetic overpressure model for 1D four-formation section with the positions and of DSA windows marked with the double arrows; b - the sensitivity chart with the major sensitivity components (HC generation constant for "Fm_2" unit) indicated with horizontal bars and their depth approximation by using the Distance Influence Function (dotted curve)

DIF(x, d) = Sj exp(-d/B,)2 (A4.2)

where xj is the j-th component of the model parameter vector x; dj is the current depth distance between the given position of depth window Z and reference i-th formation level Zt (dt = \Z - Z,|); Sj is the j-th component of the sensitivity vector S; B, - adjustable approximation parameter.

APPENDIX 5. About linearization of the effective 3D HD forward modelling operator

The combine solution of the 3D HD forward problem in basin time scale can be reduced to the consequent resolving of the following 1.5D HD forward problems for a set of 1D cross formation BM

A(z)(dP/dt) = d[C(z)(dP/dz)]/dz +8qL(z) + S(z,t) (A5.1)

and 2.5D HD forward problem solution for 2D aquifer formation BM

A(x,y)(dP/dt) = d[C(x,y)(dP/dx)]/dx+d[C(x,y)(dP/dy)]/dy+S(x,y,t) (A5.2)

with the initial and boundary conditions given in the form:

P(x, y, z, t) = 0, J d(P(x, y, z, t))/dn = 0,

d (P(x, y, z, t))/dn = P(x, y, z, t)/A j .. P(x, y, z, t) = P0(x,y z),

where P(x,y,z,t) is target over the hydrostatic pore pressure function; A(x,y,z) and C(x,y,z) are lithology controlled functions described compressibility and HD conductivity of HD rock-fluid system respectively; S(x,y,z,t) is source (additional fluid generation) term; SqL(z) - implicit Sink term given by first two right side terms of (A5.2). G = Gtu Go u Gi u Gb closed boundary of the 3D calibration area, that includes Gt - the upper surface (erosion basis) where excess pore pressure is stationary equal to zero; Gb - is the bottom part associating with non-permeable basement; Go - part of side surfaces that provides connection of the 3D area with the Remote Discharge Zone via available aquifer formations, Gi - is part of side surfaces that have not connection with the Remote Discharge Zone.

The relevant forward modelling differential operator can be defined as non-linear with respect to basis lithology constants of conduction, compaction and generation, i.e. with respect to model parameter vector x. In a final difference approximation of (A5.1-3) it is designed basing on implementation of Alternating Directions Implicit (ADI) Method (Wang, Anderson, 1982). The ADI scheme implies consequent using of four-points pattern along the rectangular 1D time-spatial grid mtz:

fflt,z = {(ti,zk) / ti+1 = ti + ht, i = 0,1,2,..n zk+1 = zk + hz, k = 0,1,2,...,nz}, (A5.4)

where ht, hz are simple time-length steps respectively.

The solutions are consequently finding for the time slice 1,2,. nt of the grid ®tz. Each particular solution is reduced to resolving of linear problem with three diagonals matrix by sweep method.

Let now the linear problem can be defined at any intermediate time slice as the following:

DiPi = di , i = 1,2,.nt, (A5.5)

where Di is square matrix of coefficients; di is vector of right side terms. The specific of three diagonals matrix D gives unique and stable ADI solution at each i-th time slice by sweep method.

The elements of matrix D and vector d are linear functions of the terms in equations (A5.1-2), which in turn are non-linear functions of model parameter x (see section 2.3). In addition di depends on Pi-1 for all

i = 1,2,.nt.

The current forward problem solution in terms of unknown over hydrostatic pore pressure vector as a function of model parameter vector can be represented at i-th time slice as the following:

Pi = (Di(x)) -1 4(x; Pi-!), i = 1,2,.nt. (A5.5*)

This recursive process along time grid should be continued till the last (present day) time slice. Thus, formally the forward modelling in terms of pressure solution can be written in form:

P = [Dnt(x)] -1 d nt(x;[Dnt-1(x)] 4 d nt-1(x,..,[D1(x)] 4 d:(x; P0).). (A5.6)

Taking into account non-linear links between basic equation terms and model parameters it is become clear that the forward modelling operator of the combine 3D HD problem is highly non-linear and non-monotonous in epy general case.

V(x,y,z)eGfe

V(x,y,z)eGfc (A53)

V(x,y,z) e Go, V(x,y,z) e Gi.