CC BY
39
The effective basin model concept and fast 3-D overpressure modelling in basin time scale
A.G. Madatov, A.-V.I. Sereda
Higher Mathematics Chair of the Polytechnical Faculty, MSTU
Abstract. The approaches and background Earth model applicable for overpressure phenomenon investigation have been reviewed and systemised in the context of their applicability for pre-drill and while-drill prediction. An Effective Basin Model (EBM) concept has been introduced basing on the successive assumption simplifying full-scale basin modelling solutions. The EBM concept allows to incorporate 1-D and 2-D approaches to overpressure modelling in basin time scale. The relevant combine 3-D overpressure modelling operator has been designed. The operator has been capable to account for the main overpressuring mechanisms including 3-D pressure communication phenomena (centroid effect, dipping aquifer beds, pressure drop across the sealing faults, etc.). A stable implicit scheme has been developed for numerical pressure solutions. The scheme reveals absolute numerical stability within the practically important range of model parameters and combined grid properties. The benefits of the developed combined 3-D HD approach versus full scale 3-D regular grid solution in computing speed is estimated to be about ^(105-106) times for the Earth model cube based on about 109 cells. The applicability of the new fast 3-D modelling scheme has been tested on real geology settings and has been checked against real data. The Gulf of Mexico basin has been used for prototyping of the centroid effect on salt induced positive Cainozoic structures. The North Sea tectonics has been used for simulating overpressured fault sealed compartments.
1. Introduction. Earth model classes in context of overpressure simulations
The wide range of Earth models were used for simulation of formation pressure and its evaluation through a long drilling history started from rather speculative (Terzaghi, Peck, 1948; Eaton, 1975) and pure empirical ("black box") methods (Huang et al., 1996; Patterson, 1996) up to highly sophisticated deterministic basin models (BM) approaches incorporating macro and local level phenomena (Ungerer, 1993; Yu, Lerch, 1996).
It is important to stress that the referred Earth models were founded to a large extent on data interpretation based on limited a priori information and the available for calibration observations usually contain their own components of errors and uncertainty. Thus, any attempt to use the calibrated Earth model for predicting purposes could give only one range of possibilities. Evidently, that the amount of source information required for calibration of the Earth dynamics system is normally invariably limited relative to the complexity of underlying excess pressure problems. It means that our pre-drill overpressure prediction is a subject to varying degrees of uncertainty. Quantification of this uncertainly is a key requirement of the drilling risk assessment and reducing.
Thus, a reasonable choice of the appropriate class of Earth models appears to be the first priority problem to be resolved on this direction.
Indeed, there is no end in producing of new and new Earth models unless their destination is determined and firm criteria of their effectiveness are introduced. The model components may be either well-founded in natural science, a crude simplification, or no more than a black box which has taken the modeller's fancy, and, equally, the prescriptions for fitting the components together may or may not be well-based. Beyond the question of the model's machinery, there are also difficult questions of the quality both of the parametrization of the model components and of the data input. Many modellers, in their enthusiasm, seem not to be fazed by these problems but push on regardless. Many appear to be possessed by their models, to have little interest in the real-world processes, and to be oblivious to the unrealities of their parametrization and input data.
In order to protect us against creating new geoscience fictions it is important to formulate the final purpose in using any specific Earth model class. Clearly, that any given Earth model class could be appropriate at the getting understanding or explanation for an investigated natural phenomenon but appears to be non effective for its simulation or prediction and opposite.
Let us schematically consider the two major classes of Earth models having potential and practically important use in context of description and prediction of overpressure phenomenon and related properties of natural geo-fluid systems. It is convenient to review them in natural sequence between pure deterministic (forward) or pure empirical (inverse) Earth models.
Producing the new Earth models for the first case can be considered as a self-sufficient research or applied task without real care about further use for data driven evaluation and/or prediction. Broadly speaking, the general destination of such Forward Models is an explanation of observable Earth phenomena and investigation of possible background mechanisms. As a matter of fact, the real data simulation is not often
considering as a primary objective of relevant forward modelling, but rather as a means of checking different hypothesis or getting deeper understanding in the relevant investigation.
As a rule, the math description of the forward Earth model includes differential equation for partial or ordinary derivatives.
In contrary, the empirical (inverse) models incorporate in our context the wide range of so-called data driven Earth models. The origin of such class of the Earth models spring from geo-engineering needs like exploration and/or production risk control, prediction, simulation, etc. Being more pragmatic by conception, the inverse models are more focused on simulation rather than explanation of observable Earth phenomena. In context of considered prediction scheme, such models include a calibration phase as a part of a model setting.
Indeed the math description of an Inverse Earth model normally implies some data fitting procedure for predefined basis fit functions. Depending on heuristic relationships chosen for data fitting this model can include elements of regression technique, data approximation approaches or highly non-linear pattern recognition algorithms.
In particular, the overpressure modelling traditionally use inverse data derived models where the main elements (normal compaction trends) were established based on regression analysis or least square data approximation (see next s/section for more details). The mimics ANN models as representatives of non-linear fitting systems are more typical for porosity and/or permeability simulation based on well log data (Patterson, 1996; Huang et al, 1996).
It is important to stress, that the essence of data driven inverse models implies the existence of sufficiently big amount of observations available for generation and calibration of heuristic relationships, whereas the shortage or even full absence of the calibration data does not really affects on effectiveness of the Forward models.
We do not eliminate a geostatistic model as a third class of the Earth model in our context. Generally, the probability theory can be applicable for both data analysis in inverse model setting and model parameter generation in forward model setting.
The scheme in Fig. 1 represents approximate position of reviewed in connection with the main topic -overpressure phenomenon Earth models on rectangular cross plot "Inverse vs. Forward models". Here, the relevant Earth models are placed in certain cells in agreement with their inherent structure and in connection with math descriptions implemented in them. Evidently, the pure deterministic forward model (Basin Models) on the one hand, and pure empirical (speculative) models, on the other hand, represent some polar-symmetric approaches to data simulation and explanation of observed geo-phenomena. The left extreme models in Fig. 1 seemed to adopt well for the prediction phase, since they include calibration phase as a part of model establishment. On the other hand, they have no room for non-formal prior geology knowledge, their parameters are poor interpretable in space and time due to speculative essence of basis data fitting functions in their background. This is why the relevant empirical relationships have quite restrict validity in the area (Dobrynin, Serebryakov, 1978; Baldwin, Butler, 1985; Campbell et al., 2000). Right extreme models could be recognised as a powerful research tool with the weak links to the real data. Thus, no one of terminal members of the Earth model family in the discussed scheme can pretend to be an optimal choice for prediction purposes.
The series of articles published by us before was aimed for development of optimally simple and suitable for data inversion purposes 1.5-D model (Madatov, Sereda, 2000a,b) and 2-D model (Madatov, Sereda, 2001) of overpressure evolution in basin time scale. Some aspects of upscaling of 3-D Earth model for the same purposes were considered in (Madatov, Sereda, 2003). In this paper we introduce a combine 3-D overpressure model as a logically next stage in development of invertable against overpressure indicators basin models. Based on reviewing of widespread background Earth models and in retrospective of achieved results we introduce concept of Effective Basin Model. This class of Earth models seems to be suitable to the relevant cell "Effective fluid dynamic models" represented in the general scheme showed in Fig. 1.
2. Review of widespread background Earth models
In agreement with the introduced above classification of Earth models implementing for abnormal pressure simulation-prediction we split appropriate Earth models on two wide categories: inverse (data driven) and forward (basin geo-fluid dynamics) models.
Fig. 1. The range of the Earth models applicable for estimation modelling of pore pressure at various scales and purposes
The difference between these two classes of Earth models is not often well defined. Indeed, some of dynamical Earth or geo-fluid models with rather sophisticated characteristics of the geological process are implementing in modern engineering approach to overpressure prediction (Swabrick et al., 1999). Naturally, they could easily be treated as a simplified BM modification (Audelt, McConnell, 1994). Conversely, upscaled basin models can apply for calibrating and can be used for empirically derived prediction. However, we believe that this subdivision is justified historically, well suits the introduced Earth model review and helps for deeper understanding the concept of effectively invertable1 basin models, which is going to be introduced and considered below.
2.1. Inverse (data driven) models
This class of the Earth models and the relevant methods of overpressure simulation-estimation-prediction historically was originated first. It incorporates empirical approaches based on explicit use of empirical relationships between the target phenomenon (excess pore pressure) and data available for calibration from surface observations (seismic velocity, interval transit time, acoustic impedance, etc.) and from wells (porosity sensitive logs, resistively logs, gas-logs, speed of penetration, etc.).
A background Earth model implemented in the relevant engineering approaches normally can be interpreted in a rather formal way ("black box" systems) or as a very gross simplification of reality. Here all internal mechanisms and causes of the formation pressure anomalies are out of consideration. The results of prediction or estimation were considered to be valuable as it is. If prediction is correct - there is no matter what a background Earth model is, if it does not work - the local empirical relationships must be rechecked against new calibration data sets.
Thus, this class of methods can formally be considered within the range of inverse (data driven) Earth
models.
A pure speculative inverse models based on identification of natural geo-fluid system with a linear (filter model) or non-linear (ANN model) casual systems as well as with pure stochastic system are not very popular in practice. A Bi-linear model for simultaneous estimation of formation pressure and lithological characteristics in interbedded sands and shales (Katz et al., 1994) could serve an example.
This model is based on series expansion of the observed well data in the form:
D(z) = Sfp(z) + SL(z) + n(z) = YB? m,
(1)
q=1
where D(z) - the well data represented in the additive model including low-frequency, the regular formation pressure component SFP(z); high-frequency, the irregular lithology component SL(z); and the random (white noise) component n(z); z - the depth; Bq - the series expansion coefficients; Wq - the suitable basis functions (usually depending on the coordinates of the layer boundary).
Decomposition of signal (formation pressure) and noise components implies least square (LS) data fitting technique for over-determined linear problem (1.1). The spatial interpretation of the set {B} inverse model parameters appears to be quite formal procedure and has had restricted validity in the Caspian Sea basin.
More complex non-linear or stochastic data fitting approaches become popular recently in connection with the reservoir simulation and history matching problems in Production Engineering (Landro, 2001, 2002). Practically, they rank intermediate positions between forward and inverse end-members in Earth model classification in Fig. 1, since they incorporate deterministic and stochastic level of model descriptions at different stage of data analysis.
The general scheme of a hybrid Earth model implemented for pressure simulation in Production Engineering can be represented as it is shown in Fig. 2.
Here, the well data is a set of well log curves and mimic system "W" can be represented, for example, by ANN permeability and/or porosity simulator (Huang et al., 1996) trained on some calibration wells. The seismic data is normally 3-D time-laps records (4-D data) available in near target reservoir locality. The mimic system "S" provides a link between seismic data and required for pressure simulation geo-fluid system (HD model) parameters: like effective stress, porosity, HC saturation, etc. The relevant inverse model can be based on Bayesian {Sen, Stoffa, 1996; Buland, Omre, 2003), System Annealing (SA) (Ingberg, 1989; Goffe et al., 1994) or Genetic (Goldberg, 1989; Mallick, 1999) seismic inversion techniques. The 2 or 3-D reservoir pressure simulation in production scale usually is based on explicit FD or FE numerical solvers of the relevant fluid-dynamic problems (Bear, Bachmat, 1991). Thus, regardless to the
Fig. 2. The global level scheme of the Hybrid Earth model implemented at pressure modelling for exploited HC reservoir (production time scale)
1 We introduce the term "INVERTABLE basin model" to emphasize the features of inverse (data driven) Earth models, which must be implemented in the forward basin model in order to adopt them, for inverting against typically available data.
Q
implemented in data inverting approach, the Earth models in the background of the relevant pressure simulators can be considered as hybrid models in agreement with the classification introduced before.
The similar hybrid models with the extent data inversion elements appear to be rather heavy from the computer point of view. Still, the inversion of windowed time-laps seismic data repeatedly available in subsurface locality of the single reservoir reveals much higher uniqueness and regularity than it does for surface pre-stack/post-stack seismic data available on whole section TWT interval.
The restrictions and stability of the surface seismic data inversion in connection with inverse Earth models are widely discussing in the special literature (see, for example, Tarantola, 1987).
Here we just conclude that implementation of a hybrid model similar to shown in Fig. 1 for overpressure prediction before drilling of exploration wells is fairly risky challenge. There are two main reasons of this: 1) the lack in formation pressure data available for calibration for whole exploring section; 2) the non-adequate data accuracy and inverse model complexity. More precisely: due to local validity of ANN methods and very unstable and non-unique solution of full waveform inversion problem at "W" and "S" system identification (see Fig. 2), respectively.
Up to now, more simple 1-D inverse models based on the mechanical compaction principles remain the broadest family of data driven Earth models applicable in formation pressure modelling-evaluation. These traditional models are generally based on empirically confirmed facts of gradual or sharp deviation of observed drilling parameters or porosity sensitive data from so-called Normal Compaction Trend lines, associated with transition from normal to abnormal formation pressure gradients (Alberty, McLean, 2003; Mouchet, Mitchell, 1989).
The typical representatives of the relevant engineering data driven approaches are the Equivalent Depth method or the Eaton method with its modifications (Magara, 1978; Eaton, 1975; Bowers, 1995; Katahara, 2003).
According to the classical Terzaghi's one-axial effective stress approximation (Terzaghi, Peck, 1948), which was empirically derived basing on massive observations of gravitational self-consolidation of marine sediments, the pore pressure P value is related to the overburden litho-static pressure Ov via the effective stress value <re as the following:
P = Ov - o-e. (2)
This ideal model ignores a tensor nature of the stress field, which includes a horizontal component. Still, according to (Stuve, 2002; Miller, 1995) this model is valid for tectonically relaxed basins, where the maximum principal load and stress are oriented vertically toward the centre of the Earth and where the same model in hydrostatic environment can be represented by the Newtonian following construction:
z
Ov(z) = Wz) [1 - #z)] - Pw(z)0(z)}dz, (2*)
0
where z - the current depth, pM(z) and pw{z) are the rock matrix and pore fluid density; <j(z) - the porosity.
It is clear that this rather general force balance equation implies instant (in a basin time scale) redistribution of the load between stressed rock grains and compressed pore fluid during the burial history. Namely, the amount of totally released pore water must gradually increase with increasing of the total sediment column weight and the relevant specific volume must be replaced with the solid fraction of the rock. This process is supposed to be controlled by the Normal Compaction Trend given empirically as a continuous function of the porosity reducing against depth (Athy, 1930) or effective stress (Smith, 1973). In case of compaction disequilibrium the pore fluid is retained and gets part of the total load on it. Consequently, the pore pressure increase increment is numerically equal to the effective stress decrease.
The described model and the relevant assumptions allowed introducing the one-by-one definitions of excess hydrostatic pressure via "equivalent depth" on porosity depth trend (Magara, 1978).
P(Z) = g[fiwZe + Pav(z - Ze)] (3)
or via Terzaghi's effective stress decrease (Eaton, 1975)
& = Onorm(^^NORM) -A. (3*)
Here, Ze denotes the shallower than z depth, at which the porosity value according to the Normal Compaction Trend is equal to the current porosity at the depth z, g - the gravity constant, pav - the average weighted rock matrix density (see formula 1.2*), a, onorm and ^/^norm - the abnormally low (due to overpressure) and normal values of the effective stress and opposite for porosity, A - the constant within the range 3-5 for loading or unloading cases, respectively.
The classical compaction disequilibrium model of Terzaghi is essentially static, i.e. does not include any parameters of geo-fluid dynamics in the basin scale. In addition, the fluid compressibility and its fraction content are ignored in it. Evidently, the validity of the relevant empirical formulas (3-3*) is ultimately limited by the case when this ideal 1-D mechanical process of overpressure generation is the paramount one among others. Whether it
is true or not depends on the relative magnitude of many factors such as a rate of sedimentation, porosity and permeability of compacting rocks, thermal environment during burial history, secondary fluid generation mechanisms, etc. (Osborne, Swabrick, 1997). Thus, the overpressure estimation based on such concept is generally suffering from the non-uniqueness problem in which too much simple trend models can fit to the available data.
The correction of the given inverse Earth model for surface erosion events implies extension of the Eaton method to additional account for the unloading compaction limb corresponded to the lower effective stress and stabilisation of the reached porosity. Practically, it means shifting and rotating of the normal trend line calibrated for the Loading limb (Bowers, 1995).
The extension for non-scalar origin of the effective stress field implies subdivision of effective (vertical) stress value in (2) on the mean stress value approximated by (Alberty, McLean, 2003):
C"mean= Ov + Ofimm + ^HmaxV3, (4)
where crmean is the mean stress; av - the vertical stress; aHmm, aHmax - the minimum and maximum horizontal stress estimations.
The vertical stress is normally easy determining from integration of density or pseudo-density logs. The minimum horizontal stress could be estimated from leak-off tests in a few in situ measurements. The maximum horizontal stress cannot be measured directly, but could be only approximately inferred from structure geology and well-bore response.
Using the mean stress instead of vertical one in Eaton like inverse Earth models appears to be obligatory within the tectonically active sub-surface environments (Yassir, Bell, 1994).
The attempts to account in the 1-D Terzaghi's consolidation model for more realistic rock rheology implies substitution of the effective stress definition (2) by the one made in the Biot theory of non-elastic deformations (Biot, 1956). The effect of inter grain plastic strain here can be accounted for by the non-linear reduction of effective stress with the depth in the form (Miller, 1995):
av = Ov - a(z)P, (5)
where a(z) is varying from 1.0 at the surface to about 0.6 at the depth 4 km2.
It is important to stress in conclusion of this sub-section, that the empirical basis for the majority of the data driven Engineering models originally were founded on the pioneer overpressure investigations in the Gulf of Mexico made by Hubbert and Rubey, 1959; Powers, 1967; Magara, 1975 and others (Magara, 1978). The specifically fast sedimentation of monotonous clay-silt sediments during the late Cainozoic obviously created here favourable conditions for application the static 1-D compaction disequilibrium model of excess pressure phenomena (Terzaghi, Peck, 1948). Quite soon, however, the classical model (2) and the relevant empirical methods required serious modifications to account for basin dynamics, heterogeneity of real section, 3-D fetures of the stress environment, etc. Still, the later application of similarly extended "porosity departure" tool for the older target intervals of the section in different sedimentary basins often reveals quite severe inconsistency. For example, in the North Sea the overpressures associated with deposition in Cretaceous-Tertiary clay rich clastic sediments do not show significant undercompaction phenomena (Thorne, Watts, 1989; Leonard, 1993; Yardley, 1998). In the offshore of Western Canada and Western Siberia significant excess pressure within the Jurassic shale-sandstone interval typically does not reveal departure from the shale normal compaction trend (Mudford, 1990; Madatov, Sereda, 2000b). Enhanced porosity in sandstones there related to the secondary diagenesis process - the dissolution of carbonate cement (Mitra, Beard, 1980).
Thus, the age of target interval of the section and the associated complexity of geology history dictate more extended empirical relationships to be established and calibrated in order to provide applicability of the data driven inverse Earth models. Clearly, that further complicating of empirical methods contradicts to their main advantage against class of forward models - a unique and quick calibration.
The list of rather popular inverse Earth models and the relevant Engineering methods with the short description represented in Table 1 gives a brief idea about wide spread concepts and approaches.
The difference between the dynamic models in Table 1 and basin models (see below) is rather conditional. Indeed, the formulas used for calibration of the relevant dynamic inverse models represent nothing else but simplified and very approximate solutions of the relevant forward problems for the geo-fluid dynamics in basin time scale given by mass conservation differential equations and equations of material state (see formulas below). These simplified formulas include some combination of basin model parameters (like sedimentation rate, temperature gradient, erosion characteristics, etc.) with cumulative properties of burying rock and released pore fluid. The framework geology model parameters normally are taken as known a priori or adjustable in the known range, whereas empirical constants require some additional calibration.
2 According to (Katsube, Caroll, 1983) a(z) = 1 - CM/C0, where CM and C0 - the skeleton and bulk rock matrix compressibility, respectively.
Table 1. The most popular engineering approaches to overpressure prediction
Reference Mechanisms accounted Short description, references Type of model
Bowers, 1995 Compaction disequilibrium affected by erosion events Extension of the Eaton method to account for possible hysteresis in normal compaction trends due to different Loading and Unloading limbs Static (basin time is not accounted) with three empirical constants
Oil and Gas Geology handbook, 1984 Excess pressure generation due to temperature dependent solid - fluid phase transformation in source rocks Fitting of the relevant kinetic model parameters to the observed gas response in drilled gas methods Oil-to-gas cracking & free gas generation according to Kinetic model
Gurevitch et al., 1987 Compaction disequilibrium affected with erosion events, thermal extension of pore fluid Extension of ED method for erosion and continuous changes of PVT fluid properties with the burial depth Static with two empirical constants
England et al., 1987 Rate of sedimentation during compaction, poor permeable vertical pressure barriers Links averaged through sedimentation column vertical overpressure gradient with averaged sedimentation rate and permeability Dynamic with one empirical constant and two geology parameters
Mann, Mackenzie, 1990 Rate of sedimentation during compaction, poor permeable vertical pressure barriers, effect of lateral fluid flow Extension of the England method to account for 3-D communication Dynamic with one empirical constant and three geology parameters
Katz et al., 1994 Undercompaction and thermal extension of pore fluid mechanisms manifested themselves in low frequency excess pressure anomaly vs. depth Join estimation of excess pressure, lithology (shale content) and noise factor based on approximation of observed pressure related data (D-exponent) by the set of orthonormal basis functions (Bi-linear model) Stochastic-Static with N unknown coefficients for the set of ortho-normal basis functions
Dobrynin, Serebryakov, 1978 Rate of sedimentation during compaction, poor permeable vertical pressure barriers, thermal extension of pore fluid Estimate abnormal pressure coefficient vs. depth as a function of sedimentation rate averaged through sedimentation column, permeability and heating rate Dynamic with four geology parameters
2.2. Forward (basin) models
The term "Basin Model" (BM) and the relevant "basin modelling" we will associate here and below with the deterministic description of geo-fluid system dynamics in basin time scale. From the mathematical point of view it is generally based on differential operators describing the global and local level processes occurring during burial and sedimentation history, i.e. in basin time scale. In contrast to this, the empirically driven methods describe the present day geo-fluid phenomena like normal compaction trends, temperature gradients, porosity-permeability relationships based on data fitting approaches. Thus, it is convenient to use this empirical basis to supply the BM determinism with the types of functional links to be implemented in analytical and/or numerical solutions and to be used for calibrations of the relevant basin models.
Originally, the BM concept is aimed for description of complex multi-stage processes of generation and secondary migration of hydrocarbons within a geo-fluid system. The relevant host sedimentary basin is described as a multi-layered porous medium affected by continuous and discontinuous deformations, heating and associated processes (Tissot, Welte, 1978; Thermal modelling..., 1985; Jouset, 1986; Ungerer et al, 1990). The temperature vs. time variation is considering here as a primarily important factor at the generation of HC. The fluid dynamic potential vs. time variation (formation pressure history) is treated as a factor which controls the secondary HC migration. Thus, rock mechanical mechanisms of overpressure generation and retaining were originally out of consideration. Still, they are implicitly included by definitions into the background deferential equations of the relevant forward models.
The basin geo-fluid dynamics grounds on the solid basis of mass and heat conservation laws and equations of material state (Lerch, 1990; Bear, Bachmat, 1991; Verweij, 1993; Miall, 2000). Depending on the generality level of the basin descriptions the relevant BM obviously involves some empirical relationships to define, for example, rock rheology, viscosity of mixed phase material, physical and chemical characteristics of deep crust and upper mantle inaccessible for observations, source rock generation constants, etc. All sub-processes involved into general geological processes at any level of their investigation reveal close mutual interconnection to each other with lots of feed back links.
The use of well grounded physical and chemical laws and assumptions indisputably create a big credit to ensure use of this approach for deeper understanding of sedimentation history and resolving present day 3-D puzzles of subsurface geological image. Still, the possibility of using basin models for prediction of the fluid system properties and, in particular, for pre-drill pore pressure prognosis with while drilling application remains a questionable problem. The reason for it is mathematical complexity of the relevant forward modelling operators and the absence of the proper link between the calibration data and tuneable model parameters. In connection with the practically important issues, it is often hard to distinguish in the BM class between a priori fixed information and the constants - terms - relationships needed to be specified additionally during calibration. Naturally, the scale of a fluid system evaluation dictates the requirements to the rank of geology setting, relationships' specification and accuracy of prediction.
Generally speaking, forward modelling in basin scale should begin from a very global level of geo-fluid system evolution, where it can be treated as a part of the developing upper crust. Being considered within the lithosphere scale, the given sedimentary basin as a whole is a current result of global level driving mechanisms creating the relevant stress and heat environments as prerequisites and controlling factors of basin tectonism and sedimentation history. Investigation of the basin model beginning from such deep roots requires application of geo-dynamics laws and plate-tectonics settings (Stuve, 2002), which are far beyond the problem of the geo-fluid system evaluation. Fortunately, the crustal stress and heat environments are to a large extent predefined a priori by the type of sedimentary basin subsidence (Allen, Allen, 1990; Miall, 2000). In particular, the Passive margins, Rift and Transform basins are created on a crust extensional basis. Extensional type basins are associated with heating of the lithosphere from below (close to vertical heat flow), stretching and thinning of the crust and close to vertical orientation of the main stress component. The sediment fill, stratigraphy patterns, compaction and sedimentation rates, faulting type and so on could be predefined in this stress and heat environment with sufficient detail to define boundary conditions for description of the relevant geo-fluid system development. From this point of view Compressional and Pull-Apart (Shearing) types of sedimentary basins seem to be less certain and suitable for the fluid dynamic modelling. The reasons of this are direction of the main stress axis which is non-vertical and
U(To,Zo)
Basin time
V
Z
Depth, Vert. Stress, Temperature
¡tí
khbbddlllll..
9*—>
t N
O
tí
№
S u(Tj,Zj)
/i
>
V
o
PH
Basin time
4 _ ^____- -
Additional fluid volume due to
- HC generation
- Thermal water extension
Pore space volume losses due to
f ™ - Gravitational compaction „.■■•** - Cementation
Fig. 3. Projection of 3-D movement of burying saturated porous rock sample U(T,D) on Depth X Time global scale pane (above) and relevant volumetric changes in pore space (draft). The loading and unloading limbs are marked with arrows, hiatus in sedimentation is marked with O mark. Note, that the hiatus in sedimentation is narrower than the one in gravitational compaction and porosity losses due to cementation are not depending on loading limb position
changeable with basin time and associated to global stress environments. The empirical links required here for description of rock compaction, heat transmission, fault architecture etc. are less certain and well proved than for the Extensional type of sedimentary basins. Formally means that it is harder to provide satisfaction of basin modelling with the necessary boundary conditions and empirical background. Still, the mass and heat conservation equations as well as fluid transport and state equations used for the fluid system analysis remain the same.
Let us briefly review the background equations commonly used at the modern basin modelling (Verweij, 1993).
These equations describe multi-mechanism physical and chemical processes occurring in representative elementary rock volume moving within global spatial-time coordinate system coupled with the developing sedimentary basin. It is important to emphasis that volumetric flows in these equations are described for the local Cartesian system coupled with the moving representative sample (see Fig. 3). Since the local system is deforming and dislocating during the
basin time scale jointly with the rock sample, the link between two systems should be introduced as well. This can be done by rewriting of continuity equations in the form of material derivatives (Bear, Bachmat, 1991) as the following:
d[M]/dt + v-V [M] = S[M]/St, (6)
where v - the vector defining current movement of the representative sample of rock coupled with local coordinate system relative to global basin system; [M] could be the mass or volume of the representative sample. For example, the conservation of solid mass during burial history can be expressed as follows in the global (stationary) system:
8W\ - 0]/St = 0. (6*)
The same in the local (movable) coordinate system is given by:
d[p0(1 - $)]/dt = -Vq0[^0(1 - 0]. (6**)
The vector v controlling the local system movement is a commutative response on global scale dislocations of representative volume involved into macro-level processes like: crust subsidence, sediment deposition, erosion, tectonics, etc. Normally it is represented by a priori geology input and accounted in the form of boundary conditions, assumptions and process specifications according to the predefined type of sedimentary basin, its stratigraphy pattern, tectonic signature and so on. Moreover, it could be stated that knowledge and assumptions about this vector define to a certain extent the class of the basin model. Globally controlled conservation of mass:
{ d[p0(1 - <P)]/dt = V-[00(1 - 0q0] - Q (7)
dfafi/dt = V-foqO + Q
where q0 and q1 are the solid (0) and fluid (1) phase volumetric flow speed relative to the global coordinate system; the specific volume (porosity); ,00,1 - the solid (0) and fluid (1) phase density; SQ0,i - the solid (0) and fluid (1) source terms (SQ0 = <5Q1).
This mass conservation equation system describes in a rather general form (Shi, Wang, 1986) balancing between the solid rock and single-phase composite fluid mass inside a virtually uniform sample volume U during its burial history. The densities of the solid and fluid parts are treated as global time and space depending continuous functions.
The implementation of Terzaghi's model (2) and lineariation of time derivatives for specific volume and phase densities allows reducing (7) to the single differential equation (Madatov, Sereda, 2000a):
[C/(1 - <jj) - <p(c1- c0)]dP/dt = C/(1 - <f>)dOv/dt + [B/(1 - $ + <Kfi1-/30)]dT/dt + SG1 + div(q0, (7*)
where C, B - the bulk compressibility and expansibility of the porous media; c0,i - the compressibility of solid and fluid components, respectively; - the thermal expansibility of solid and fluid components, respectively; SG1 -the specific storage of additional pore fluid volume due to the difference between the densities of converting solid to fluid phases; P and Ov - given by (2); dT/dt - the temperature changes controlled by the heat transfer process (9). At that, any phase density is expressible in the linearised form (Domenico, Palciauscas, 1979):
(1/p)(dp/dt) = c (dP/dt) + ¡3 (dT/dt). (7**)
The rock matrix deformation is treated as irreversible porosity losses defining in agreement to the effective stress increment - Smiths' compaction model (Smith, 1973) and temperature increment - the secondary porosity losses model (Rittenhouse, 1971)
(1/$(d^/dt) = C (da/dt) + B (dT/dt). (7***)
Locally controlled conservation of solute fluid mass:
- V-r,- V-{(q/$ s,} + Ssj = dsj / dt, (8)
where rj is the speed of specific mass changes for the j-th solute component transmitted through a single area element orthogonal to the given component of the concentration gradient; sj - the specific mass of a single chemical solute per unit volume of the aqueous solution (instant concentration); Ss, - the relevant net source (+) or sink (-) term.
The mass conservation described by (8) implies balancing mineral solution and phase transform against mineral precipitation and internal degradation controlled by the local space scale variation in the component concentration given in the basin time scale.
Conservation of heat: J J J
V[@0 (1 - <P) + <PT Sj0j- £ j Sj ]VT + ST = [(1 - <PT SjYj] dT/dt, (9)
j= 1 j=1 j=1
where 0- the thermal conductivity of the solid (0) and the j-th components of the pore fluid; Sj - the fraction of the j-th saturated fluid component (saturation); ST - the heat source (+) or sink (-) term; y - the heat capacity of the solid (0) and the j-th components of the pore fluid.
Flow steady state transport equation (generalized Darcy law):
q- q0= DPj^P + Dt,VT. (10)
where qj is the steady state speed (volumetric flow rate) of the fluid j-th phase relative to the solid phase volumetric speed q0; Dp, and DTj are the pressure and temperature (phenomenological) induced conduction coefficients. At that, the conductivity ofj-th fluid components is determined as the following:
DPj = K^ /#, (11)
where K0 - the rock matrix permeability that could be a tensor in the general case, Kj - the relative j-th phase permeability (scalar); ^ - the viscosity of the j-th components of the pore fluid.
Diffusion steady state transport equation (Fick law):
r= FVsp (12)
where r is the steady state speed (rate) of specific mass changes in agreement with (8) definition; F - the coefficient of diffusion (can be a tensor); s, - the mass of a single chemical solute per unit volume of the aqueous solution.
Equations of material state for fluid and solid phases:
PVT functions define the density and viscosity of each j-th component of the pore fluid as smoothed functions of current pore pressure, temperature and fraction content (Denesh, 1998):
{ p3 = flP,T,tSj); (13)
(j = fi(P,T,$Sj).
Empirical trend functions define combined mechanical (rheology) and chemical (diagenesis) responses of the sediment matrix (solid skeleton) in terms of the representative specific volume (porosity) and pore space geometry (permeability) characteristics. The porosity and permeability as an effective measure on matrix compaction and hydraulic conduction of rock imply substitution of real heterogeneous porous media on pseudo-homogeneous solid material with the empirically derived rheology.
{ ^ = /Cc,cr,T,Z); (13*)
* K, = rts,0, 1 )
where c, s - the rock matrix compaction and specific surface constants given in the vector form.
The implementation of this theoretical background accommodated to the boundary conditions follows from the basin type and available empirical relationships allowed to create the wide range of deterministic quantitative models (forward Earth model) describing fluid dynamics in basin time scale (Allen, Allen, 1990; Tissot, Welte, 1978; Guidish et al, 1985; Doligez et al., 1985; Yu, Lerch, 1996). The relevant set of a final difference or final element operators allow simulation of: burial and thermal history; sediment compaction; retention and generation of excess pore pressure; generation and migration of hydrocarbons (HC) in the form of single- and multi-phase fluid flows; trapping, accumulation and integrity of possible HC plays and so on. An extended and exhaustive review of available basin modelling approaches and packages is available in the published literature (see, for example, Allen, Allen, 1990; Lerch, 1990; Miall, 2000).
Most of the modern multidisciplinary BM tools are mainly focused on primary and secondary migration of HC during developing of geo-fluidal Earth system. This complex physical-chemical process incorporates a set of global level processes occurring in nature simultaneously in terms of the basin time scale (Guidish et al., 1985). This set incorporates the following disciplines and related processes:
• The geo-dynamics in the part dealing with the upper crust tectonics, rock material rheology, flexure and disjunctive deformation in connection with structure evolution compensating lithosphere blocks' movement.
• The sedimentology which describes the burial history of sediments as a consequence of processes to include: sediment deposition, accumulation, and consolidation in connection with structural evolution compensating basement subsidence.
• The thermodynamic of the global and local sources and heat transmission through the non-homogeneously conductive porous media as a thermal history of sediments.
• Organic geochemistry controlling generation of HC from source rock - primary oil and gas migration.
• Multi-phase underground fluid dynamics driving upward and lateral flow of dissolved and free HC components through the porous media - secondary oil and gas migration.
I
Fracturing
Stress
y
Burial/ Subsidence and Tec tonics
Porosity and permeability
Fracture filliij»
y permeability ¿
Chemical compaction and diagenesis
Pore Du id pressure
M el han ugenesis
X
Fig. 4. The feedback loops suggesting the many coupled processes operating in sedimentary basin (according to RTM Basin Simulators, Ortoleva, 1999)
Clearly, that each of these global level processes include implicitly a set of the underlying parallel local level processes, where each of them can also be represented by the secondary local level process combination and so on. For example, the sediment accumulation includes basement motions, eustatic sea level variation, compaction of sediment matrix, release of pore fluid, faulting erosion, etc. The compaction itself can represent as a series of simultaneous physical and chemical processes of loosing or conservation of pore space in burying rocks (Waples, Kamata, 1993).
Finally, the attempt to account for at least primarily important basin processes leads to highly non-linear forward model with the multiple enclosed feedback links (see Fig. 4).
The applicability of this kind of basin simulators even for rank prediction of HC accumulation, present day integrity, content, size, etc. is problematic.
Indeed, the presently available data input and geology knowledge are results of data transformation, collection and interpretation of a priori information which in the best case can be treated as a present day imprint of a multi-level process. The geology image as an external framework is certainly formed by a host of different processes that have occurred in an infinite variety of combinations throughout geological history. Most of the regional scale geological conditions that we deal with consist of strata that have been initially formed with relatively simple geometry by depositional processes. The shape complications have then been introduced by later continuous and discontinuous deformations, erosions, intrusions and other tectonic activity. Strictly speaking, any 3-D geometry reconstruction back in basin time requires some paleo-markers to be available in order to follow unique way. Normally they are implicitly substituted with intrinsic assumptions or hypothesis. Even in case when such paleo-markers can be made available from the present day observations, like palaeontology indicators (bio-markers) or vitrinite reflectance, the relevant deposition environment and reached underground temperature, respectively, is a matter of interpreting and object of errors.
Today it is commonly accepted that an attempt to use such kind of deterministic process simulators for prediction purposes can cause serious limitation on the mathematical complexity of the corresponding forward model (Lerch, 1991; Zhao, Lerch, 1993; Yu et al., 1995; Madatov, Sereda, 2003).
What is the optimal complexity of a basin model implemented for overpressure prediction? What is the strategy of its simplification?
It is intuitively clear that the way of simplification consists in optimal choice of workable dimension of basis continuity equations (6-9), introducing of assumptions that limit intrinsic ambiguity of multi-mechanism processes and links with the data and a priori information capable to constrain model uncertainty inside of the chosen class of models. The only formal way of doing this is formulation of the relevant calibration task for the forward BM in terms of the relevant Inverse Problem solution.
With respect to prediction of HC generation, migration and accumulation the Inverse Problem for dynamical indicators of burial and thermal history was first formulated by Lerch in 1991. The principles of inverting of excess hydrostatic pressure indicators were first introduced by us in 1994 and then developed in series of publications (Madatov et al, 1996; 1997; 1999; Madatov, Sereda, 1997; 2000a). On a qualitative level they were aimed to coordinate a particular forward model concept (underlying assumptions) with a priori available geology information (type of sedimentary basin, local stratigraphy-lithology tectonic settings, etc.) (Madatov, Sereda, 2003). On a quantitative level they determine the type and dimension of the relevant forward operator, specific model space for inverting and top criteria of inversion routine in agreement with the calibration data quality (Madatov, Sereda, 2000a).
Now we will focus on effective basin model concept in connection with aimed applicability of the relevant forward modelling operator for inversion with respect to available overpressure indicators.
3. An Effective Basin Model concept in context of excess pressure prediction
The concept of Effective Earth Models (EEM) is widespread in geophysical data interpreting and inversion (Berdnikov et al., 1985; Claerbout, 1985; Helth et al., 1994). In seismic data inversion it implies substitution of the inherently inhomogeneous media by the apparently homogeneous prototype on "INPUT -OUTPUT" level. The level of response depends on the degree of uniformity acceptable for observation and simulation of a target geo-phenomenon.
Even the most uniform material is made up of atoms, so its properties are grossly non-uniform when viewed on a small enough scale. The material made of obvious structural elements may be highly uniform when viewed on a large scale. If choice of scale were entirely arbitrary, the term homogeneity would not be particularly useful. The context always provides some reference uniformity level to serve as a scale and accuracy of measurement and prediction. From the System Analysis point of view that means replacing of multi-level
highly non-linear links established between parameters in some basis deterministic model with the set of few linearly independent ones without losing the main functional features addressed by the investigation. For example, the Effective Earth Model of the same target Earth volume could be represented in rather different way on the seismic data inversion depending on prediction purposes. When the target is a Structure image, the best EEM is defined as a stack of 3-D reflector sub-surfaces. When the target is a local acoustic Impedance image the best EEM is described as a 1-D high frequency reflectivity response (ARMA time series).
Regardless of the purposes of prediction and forward model of phenomena the principles of effectiveness of model parameters remain the same:
1) Any internal and external functional links established for the coefficients determined in the basis equation(s) within the given class of the forward models can be effectively replaced by constants or more simple links with fewer constants to be defined if the relevant forward response of the model can simulate all possible data vagaries with the predefined accuracy. This lowest-useful set of coefficients to be specified for simulation and calibrated for prediction defines a Control Model Parameter set in effective representation.
2) Each of the Control Model Parameters must have cumulative representability, i.e. it must reveal as uniform as possible the behaviour inside the given intrinsically homogeneous Earth model unit and be distinguishable against an adjacent unit parameter.
3) The Control Model Parameters must be interpretable in the context of assumptions introduced for the given forward model class.
Let us now review the basin model approach in order to adjust it for overpressure prediction purposes, to introduce necessary assumptions and to specify the sufficiently complete Control Model Parameter set for the Effective Basin Model.
Before reviewing it is important to emphasis some general requirements to the Effective Earth model to be satisfied to the discussed topic - a 3-D overpressure simulation useful for prediction purposes. Let us formulate them as the following:
1) The EEM should describe a target phenomenon in the simplest and unique way.
2) The EEM must have places to include as much a priori knowledge about the target phenomenon and related to it data sets as possible.
3) The amount and specification of data required for calibrating of EEM must be well balanced with the real capacity that you could afford at the data inverting.
4) The data related to the predicted phenomena have to be linked to the EEM in the simplest and unique way.
5) The prediction routine must be continuously sensitive to the amount and quality of available for calibration data. In other words, the uncertainty (risk) assessment shall be part of prediction and shall converge to minimum at increasing data involved into calibration.
6) The EEM and prediction routine must support real time updating as more calibration data is made available during drilling.
It is clear that aimed Effective Basin Model (EBM) reveals the features of both inverse and forward Earth models considered above.
The benefits of using the deterministic forward model in the BM approach as a basis for data inverting consists in the following possibilities:
• including a priori knowledge about local geology settings in compact and deterministic way;
• explicit account for different hypothesis about overpressuring mechanisms;
• increase in understanding of origins and magnitudes of presently acting overpressure phenomena.
The disadvantage of using the BM approach for pore pressure prediction follows from the fact that BM tools were not created for prediction purposes, but for modelling. In particular, the primary objectives of the BM approach are simulation of HC primary and secondary migration in basin time retrospective. Thus, the Earth model must first be upscaled for pore pressure purposes, i.e. be prepared for inverting available pressure. The approach to the Earth model upscaling in the discussed context was considered in details in (Madatov, Sereda, 2003).
Our task now is more general. We have to start from a very general BM level given in (6-13*) and consequently simplify it keeping balance between reasonable complexity of a forward overpressure model in basin scale and practically acceptable accuracy of overpressure prediction.
The overview of the main components of basin models will help us to specify deterministic background for the aimed effective forward overpressure model and start its construction based on consequent BM simplifications. The first step in such simplification is recognition of the main well proved overpressuring mechanisms among the set of simultaneously occurring the global level processes. Then they have to be ranged in agreement with sensitive overpressure indicators and linked via the model parameters with the corresponding terms in the basis of the basin modelling equations (6-13*). The next step is subdivision of the coupled processes into primary and secondary important and introducing the relevant assumptions and limitations to allocate the suitable classes of effective basin models.
To begin with, it is convenient to represent the main basin model components in speculative deductive sequence: Processes > Responses > Phenomena as it is shown in Fig. 5. Here, each link arrow is colour coded3 according to the level of importance of the relevant influence from red (most important) down to blue (least important). Note that in the nature each component of the determined classes often has casually mixed with adjacent class origin and firmly linked with them by feedback relationships like it is shown in Fig. 4. Indeed, fast subsidence and sedimentation create condition for the fast heating of sediments, but accompanying this process compaction disequilibrium preserves the abnormally high porosity level, that leads to significant decrease of thermal conduction. Finally, the quickly buried areas associated with abnormal cooling zone in 3-D heat flux field (the blanketing effect) which is impossible to explain without the mentioned feedback links. Thus, some of the arrow links could have an inverse arrow mark as well.
Note, that subdivision of Processes and Responses on three and four separated groups respectively is conventional speculative construction made in order to follow common viewing accepted in the BM approach and designing of the relevant BM tools. In particular, the burial history is inseparable from tectonic and thermal environment. Moreover, in nature it controls them and is controlled by them at the same time. Here, we split subvertical (across the beds) and sub-horizontal (lateral) stress environment during basin evolution in order to follow common basin types and explicitly specify them via different terms in basis equations. Thermal history is treated here as a separate part of the sedimentation history to highlight corresponding links with mostly temperature dependent secondary processes (Responses) and the overpressure mechanisms induced by them.
Let us now briefly consider each of the processes in the context of basic BM equations (7-13). It will help us to allocate their relative contributions in the excess pressure generation and/or retention.
3.1. Tectonic compression
The deviation from Extensional Type of sedimentary Basin during evolution of the relevant geo-fluid system as a whole or local lateral compression in period of fault or intrusion activation can create essentially non-vertical confined or local stress environment, respectively.
Simulation of lateral stress field dynamics at basin scale on global and local levels is matter of Geodynamics of the lithosphere plate and single basin level, respectively (Stuve, 2002). Formally, this process is represented implicitly by the vector v in (6) controlling the local system movement in the global time-spatial coordinate system. The lateral components of the current stress increment in basin time are expressed by dajdt + dajdt in mass conservation equations (7-7***).
The tectonically induced stress term in the global scale mass balance equations (7-7***) can be treated effectively as an additional to daz /dt component increasing rock matrix compaction. Based on this viewing, it will cause the overpressuring within some zones where the corresponding additional volume of pore fluid cannot be released during fixed time of burial history. If this fixed time gate is limited by present day, then excess pore pressure will be a matter of certain drilling risk. A good example of this phenomenon gives a wide overpressured zone in front of San Andreas fault in California (Morley, 1999). Overpressure build up due to this origin can be very rapid. Thus, it can be especially important for the basins that are tectonically active up to the present time. The problem of estimation and quantitative prediction of the overpressuring related to this process is difficult due to poor understanding of all underlying global processes and complexity of the relevant geology settings (Osborne, Swabric, 1997). Locally, it can create prerequisite for lateral fluid barriers and overpressure compartmentalisation due to generation of the sealing fault systems. Each single fault with non-zero displacement and low membrane transmissibility can act as both pressure drain (at the very beginning of strata displacement) and pressure seal (in later time) for a laterally permeable aquifer formation (Antonellini, Aydin, 1995; Knott, 1993). The overpressure mechanisms induced by this local tectonic activity can be recognised as an essentially secondary process of 3-D excess pressure redistribution.
3.2. Burial history
In this context we define burial history as a part of the global level processes responsible for motion of basement - sediment fill, paleo-bathymetry and eustatic sea level variation vs. basin time. This limitation
3 A colored version of this article is available in electronic form in web-site of the given issue (http://vestnik.mstu.edu.ru)
Fig. 5. A deductive sequence: Processes > Responses > Phenomena reduced to Overpressuring phenomenon with the arrow color indicating sensitivity of link (1 - low; 2 - middle; 3 - high)
reduces burial history to 1-D model of basin bottom subsidence - uplift and the relevant deposition - erosion of sediment column against depth axis. Such model appears to be well developed and supplied with proven empirical relationships (Allen, Allen, 1990; Miall, 2000).
The robust model of compensated deposition implies that the sediment fill is balancing by basement subsidence. At that the dipping of the basin bottom is compensated by further sediment deposition and sea level plus water depth variations. Conversely, the uplift of the basin bottom is accompanied by some erosion or depositional hiatus (Guidish et al, 1985; Barker, 1996). Thus, the vertical loading-unloading process and corresponding stress component daz /dt in the mass conservation equation (7*) is controlled by vertical basement motions. There are at least three causes for movement of the basin bottom: isostatic subsidence, flexure subsidence and thermal subsidence (Miall, 2000; Stuve, 2002). Which one of the relevant deterministic geo-dynamic model has to be implemented? How to determine the required coefficients in basis equations in order to describe the corresponding lithosphere compensation model? These problems again create an uncertainty in defining the burial history by using the forward geo-dynamic problem approach.
The 1-D (along depth axis) treatment of burial history accepted in our consideration (Madatov, Sereda, 2000a) allows implementation of the inverse 1-D deposition model and corresponding backstripping technique (Allen, Allen, 1990) to reconstruct the vertical component of basement movement at any given X, Y location within the basin (see Fig. 6).
^¡iS? He
a b c
Fig. 6. The three types of stress environments for different basin genesis: a - diverging basins; b - converging basins; c - shearing basins (after Miall, 2000)
The inverse approach at the burial backstripping is based on assumption about irreversible porosity reduction during rock matrix compaction and conservation of the solid part of sediment during continuous burial (6). It delivers the unique and certain decompaction operator under the conditions that compaction trend for the given lithology can be calibrated against available porosity data4 (Schneider et al., 1994; Liu, Roaldset, 1994) in case of continuous loading of the sediments. Erosion events in 1-D backstripping approach are uniquely linked with upward basement motions regardless of their lithosphere origins. Still recovering of basement uplift by using the inverse 1-D erosion model remains less certain problem.
According to I. Lerch (1990) there are three possibilities to build the inverse 1-D erosion model based on time interval of the detected unconformity: (1) as a depositional hiatus; (2) based on extrapolation of sedimentation rate forward from below formation as an loading rate and back from above formation as an unloading rate on condition that collected and eroded during given time gape thickness is the same; (3) based on local geology knowledge. Additional information can be extracted from so-called maximum reached depth markers. These paleo load and/or temperature markers are linked with porosity vs. depth data collected for the given formation unit. Indeed, porosity losses established for given lithology according to its loading limb deliver one-by-one link with maximum reached depth/effective stress on condition of irreversible behaviour of compaction (Reike, Chilingarian, 1974). The quantitative methods based on 1-D well data analysis (Reike, Chilingarian, 1974) and 2-D seismic-stratigraphy (Seismic Stratigraphy, 1977) are well documented.
4 TThe initial thickness of every lithologically uniformed formation unit can be recovered (inverted back to the start of its burial history) uniquely basing on the present day thickness (hT) and the averaged porosity (<jr ) observations and the initial porosity value (^0) available from the trend as following: hTh0 = hT (1-^/(1-$)).
p
200 175 150 125 10Û 75 50 25 0 Time
Fig. 7. The 1-D Inverse (back stripping) formation model of burial history for a single well with subsidence (pink) paleo-bathymetric (blue) and eustatic sea level curves attached
The burial history in the given 1-D stress environment context controls mainly compaction disequilibrium as a mostly well proven origin of overpressuring. From this point of view the relevant backstripping technique can be considered as an extension of the classical Equivalent Depth or Eaton static models (see s/section 2.1) used up to date in the engineering approach to overpressure prediction. Overpressure transmission and sealing mechanisms can also be consequences of the differential rate in vertical stress development during 1-D burial history. Indeed, a significantly different burial rate for the same aquifer formation over a short X, Y distance can result in lateral pressure transmission and appearance of Centroid-like effect in the arc part of positive structures (Traugott, 1996) (see also the modelling examples below). The vertical compaction of poor permeable shale induced by rapid sedimentation rate and loading can lead to creation of vertical pressure barriers (seals). Jointly with increasing the underground temperature, this process induces mineral transformation like Smectite to Illite or Gypsum-to-Anbydrite that in turn can result with the local overpressuring due to dehydration and reducing permeability (Osborne, Swarbrick, 1997).
3.3. Thermal history
The model of temperature evolution on a basin scale normally implies introducing of two global level processes: heat generation and heat transfer. Their balance on the basin time scale is described by the relevant heat conservation differential equation (9). The heat generators as a thermal input can be treated, similarly to the tectonics input, on the lithosphere (supra-basin) level and sedimentary cover (single basin) levels. There are several basal heat flow models describing heat generation on the global level: McKenzie's lithosphere stretching model, Sleep's thermal expansion model, Falvey's deep crustal metamorphism model and so on (Stuve, 2002). Some of them like Royden's dike intrusion model have universal applicability either on global or local levels (see for more details Thermal modelling in sedimentary basins, 1985).
It is commonly accepted that the global level basal heat flow has two components: radiogenic contribution of the crust and sub-crustal heat flow (Ungerer et al, 1990). The sub-crustal heat flow is simulated as a constant (Foreland basin type) or smoothly decreasing (Rift basin type) in post activation time interval. The radiogenic part is estimated as an addition to sub-crustal heat flow generated within the radioactive crust thickness according to Radiogenic Production function empirically defined basing on seismic velocity analysis (Rybach, Bunterbarth, 1984).
The heat conduction and in less degree the convection are commonly recognised as dominant processes of thermal transfer in sedimentary rocks. The thermal conductivity of rock controls the slow and smoothed variation of the temperature gradient in depth and in basin time. This property of the rock is highly sensitive to porosity and is normally defined from empirical trends specified for particular lithology (Chapman et al, 1984). Typically two or three empirically based constants are required to be calibrated against paleo-temperature markers and present day temperature data from wells and from surface observation. The predominant direction of the heat transfer is subvertical. However, some lateral components could be related to the flanks of rift basins on the global level, magmatic intrusion, salt domes and crystalline block to cause lateral heat transfer on a local level.
The temperature dependent overpressuring mechanisms can be classified as primary and secondary. The main primary mechanism is known as aqua-thermal pressuring (Shi, Wang, 1986). It is caused by volumetric expansion of pore fluid with increasing temperature and is controlled by its phase content and PVT properties of each fraction according to current saturation (8) and equation of state (13). This effect is very small in comparison with volume losses during compaction. For example, fresh water heated from 50°C up to 95°C has a relative increase in volume of less than 2 %. Thus, the aquathermal pressuring effect can generate significant overpressuring effect only in a completely sealed condition. However, such conditions are very unstable in situ due to fracturing effects, occurring immediately after approaching of inter-grain stress to the relevant micro-cracking limit (increasing excess pore pressure up to fracturing limit). Luo and Vasseur (1992) showed that the contribution of this primary temperature effect is negligible down to depth 5-6 km at an average temperature gradient below 35°C/km.
The overpressuring mechanisms related to the diagenetic processes and HC generation are treated here as secondary temperature-dependent and considered in connection with the relevant Responses (see Fig. 5). Despite of the secondary level these mechanisms can play a paramount role in local interval overpressuring. However, the effects of additional pore fluid generation on exceeding hydrostatic pressure regime is less stable in basin time scale, sharply changeable with local source rock in situ conditions, porosity and permeability environments and so on. The related forward modelling and prediction problems are also discussing in s/section 3.5.
3.4. Mechanical rock response
The mechanical response of the sedimentary rock matrix continuously stressed by overloaded material and compressed by additional tectonics components of confined stress is reduced to plastic and fragile deformations on a micro level (Charlez, 1991; Straughan et al, 2001) or to flexure and disjunctive dislocations on a macro level (Miall, 2000).
The micro-level response can be understood as a local accommodation of rock to plastic deformation. Starting from elastic deformation it consequently goes through grain rearrangement and pore crushing up to viscous flow in response to gradual increase in confined stress level achieved by sample volume during burial history (see Fig. 3). In the context of consideration and creation of a forward model of the process, it is important to define a strain parameter, stress factor and rheology model. Generally, the BM approach (for example, in the TEMISPACK model (Doligez et al, 1985; Ungerer et al., 1990)) uses: porosity losses as an effective parameter of strain for pseudo-homogeneous porous media, Terzaghi's one-axial effective stress model (2) and rheology based on the model of irreversible porosity losses following a loading limb. Clearly, such a concept implies that horizontal components of confined stress are negligible and that the small volume of representative sample remains preserved during burial history. This model could be recognised in terms of gravitational compaction or consolidation (Magara, 1978). In connection with the pore fluid saturated rock sample in Fig.3 it is associated with its irreversible expulsion from continuously decreasing pore volume.
There are various consolidation models developed in soil mechanics, rock mechanics and sedimentology (Reike, Chilingarian, 1974). In connection with the pore pressure prediction approaches the classical 1-D model of saturated clay compaction (Terzaghi, Peck, 1948) is still used as a basis in many Engineering approaches (see Table 1). A typical example of a highly effective model that ignores the underlying consolidation processes and mechanisms in return for getting fast and rather certain output in target scale is the Smith (1973) formula:
<j) = mj[exp(- m2cr)], (14)
where consolidation constants m1 (initial porosity) and m2 (compaction constant) are objects of calibration against borehole data.
On closer investigation, mechanical compaction includes at least three simultaneous micro-scale mechanisms: re-packing, crushing and welding (Schneider et al., 1994; Waples, Kamata, 1993). The pressure solution can be considered as adjacent to chemical response (see below).
The most general consolidation process is repacking, which implies rearrangement of solid particles (grains) during increasing confining load from initial to more compacted structure. Repacking is expected to begin from about 100 m burial depth and to decay by about 2500-3000 m (Moore, 1989). The empirical relationships suggested for description of repacking required typically two compaction constants to be calibrated: in exponential (Athy-like formula) or power law formula (Baldwin, Butler, 1985).
Crushing of weak and hollow grains' mechanisms are practically important for carbonates (especially for chalks), where they generate so-called porosity collapse at certain stress level (Reike, Chilingarian, 1974). The relevant empirical relationships (Moore, 1989) are available for the final difference porosity increment form of the exponential trend function. In addition to two trend compaction constants they require one more collapsing constant to be calibrated.
The reduction of porosity due to ductile flow is a rheology process common for carbonates and sandstones, and less common for shale. The empirical relationship for porosity loss is given by (Mitra, Beard, 1980) in the recursive form:
$ = 1 - Vr/( Vr + Vp + A Vd), (15)
where Vr, Vp - the volume of solid and pore spaces, respectively, prior to ductile flow; AVD = - CD (D - 0.02AZ), CD - the constant; D - the fraction of ductile grains which do not yet suffer from ductile flow (empirical constant); AZ - the depth increment.
Thus, a process-oriented approach to evaluation of porosity losses due to mechanical compaction involves at least three underlying processes and requires at least six empirical consolidation constants to be calibrated.
At the same time there is only one coefficient in the basis mass conservation equation (7) to control volumetric changes between solid and fluid fraction within the representative rock sample. It is bulk compressibility of the porous media c0. In agreement with (7*) this coefficient can be considered as a cumulative Effective Basin Model (EBM) parameter uniquely specified for the given formation. Evidently, there is an infinite number of options for the functional representation of rock bulk compressibility as an effective controlling factor of porosity reduction: from a single end-member constant for hole section up to a function of all underlying compaction sub-processes with dozens of end-member empirical constants for every terminal sub-process.
The macro (tectonic) level response can be treated in terms of basin structure evolution: folding, clay and salt diapirism, etc. Except for exotic cases like a wrapped dome, this response is associated with the activity of lateral stress components. The evolution of the relevant basin level models was briefly discussed in s/section 3.1.
Let us now briefly consider the mechanical response of sediments which exceed the plastic deformation level. The relevant disjunctive dislocations on macro-level and fragile deformations on micro-level imprint themselves in the present day fault systems (Fault Mechanics..., 1992) and fracture network (Domenico, Palciauskas, 1979), respectively.
The typical micro-level fragile response of rock matrix is due to pore pressure exceeding fracture limit. In both cases fracturing is an over-plastic response of deformed rock material on strain. The occurrence of faulting and fracturing can be predefined based on paleo-reconstruction of the stress environment in agreement with the rock mechanic principles (Straughan et al, 2001). One of the common approaches consists in defining the "ductility" as an upper limit of the strain to be accommodated by plastic deformation. Reconstructing or forward modelling of total stress environment on global and local levels allows detection of time and place when and where the strain event exceeded this threshold and faulting or fracturing was activated respectively. The popular Normal Ductility trend formula is given in the form (Charlez, 1991):
DUCTILITY = Klexp(K2STRESS). (16)
This relationship can be established and calibrated for each specific lithology as a smoothed function of depth by analogy with the Normal Compaction trend. The number of trend constants (Ki, K2) required to be determined/calibrated depends on the detail of the stress field model. As a minimum this link must be established for compressive and extension strains like it does in the Virtual Basement Displacement (VBD) method, or in the Basin2 method (Rowan, 1993). Much more parameters are required to be entered in RTM FRAC BM simulator (Ortovela, 1999) that implements, for instance, tensor stress versus pore pressure fields to simulate growth, geometry orientation and intensity of complex fault and fracture network on macro and local levels of the forward modelling. These deterministic essentially dynamic models of the over-plastic rock response allow investigation of the pulse of fault activity, cycles of pressure release via opened micro-fractures and pressure seal caused by healing fracture network. Clearly, neither calibration nor even verification of such complex deterministic models with extended input requirements are practically solvable problems. By analogy with flexure deformation models, the empirically derived models with much fewer constants to be adjusted/tuned appear to be robust for evaluation of fault transmissibility and fracture permeability as end-members of the effective phenomenological models on macro and micro levels, respectively.
Both plastic and fragile responses of rock matrix influence on rock transport properties on the local scale level. These transport functions of the rock mechanical responses are represented in Transport (10) and Heat Transfer (9) equations by effective coefficients: DPJ and DTj - the fluid conductivity and 0- the thermal conductivity. Crucially important in pore pressure prediction context is rock matrix permeability K0, which is a combine function of mechanical and chemical responses on pore space fabric, micro-channel connectivity, tortuosity, etc. Normally, pore permeability as a function of rock flexure mechanical and chemical responses and fracture permeability as a function of rock disjunctive mechanical response are distinguished.
The pore permeability is normally determined as a smoothed function of porosity (Reike, Chilingarian, 1974; Magara, 1978). The type of function link and range of necessary constants are predefined empirically for each given lithology. There are several well proven relationships established for clastic rocks and carbonates. For example, Koseny-Carman relationships for mudrocks (the North Sea):
Ka = 0.2 ^3/[S2(1 - #], (17)
Mann and McKenzie (1990) for sandstone (the Gulf of Mexico):
Kss = S"V, (17*)
Schole (1977) for chalks (the North Sea):
KCh= 10(6.25^-4.68) S"2, (17**)
where K - the permeability [m2]5 with lower index for lithology; <j) - the porosity (in fractions); S - the specific surface constant in m-1 that defines internal pore space fabric for the given lithology.
There are also several robust models of the Fracturing permeability based on similarly simple empirical relationships. For example, in the (Rowan, 1993) model the effective fracture permeability is defined as a function of the stain threshold SD:
Log(Kef) = K3/ Sd + K4, (18)
where SD derived from the Normal Ductility trend for specific lithology includes two empirical constants in addition to K1, K2 introduced in (13*), K3, K4 - the empirical constants.
The fluid transport and heat transfer properties on the global scale level are controlled by faulting. Relevant fluid and/or heat transmissibility of the fault zone are part of the global scale environment, that indirectly affect the rock sample and relevant local scale pore pressure model (see Fig. 3). These elements of the framework basin model could be evaluated as a part by product of the disjunctive mechanical response simulation on the basin scale. However, this way is rather non-practical due to the problems with verification / calibrations of the relevant complex deterministic forward models of tectonic evolution mentioned above.
5 1 Darcy = 10"12 m2.
Most of the BM tools implement the inverse model to account for the fault transmissibility, i.e. allow the user to recognise the fault on sub-surface structure map and introduce it rather then predict it appearing and evolution based on geo-dynamic principles. This essentially effective approach assumes simple predefined dynamic of the faulting process and allows calibration of the control hydraulic parameters based on the phenomenological model of the fault transmissibility. For example, the following empirical relationships suggested by (Manzocchi et al., 1998) are used to define combine fault permeability as a function of the shale content and fault displacement:
Log (Kf) = - 4 SGR - 0.25 Log(D)(1 - SGR)5, (18*)
where SGR - the representative measure of the shale content (Shale Gouge Ratio) given as a fraction. The simplified fault activation model treats faulting as a catastrophic event (Madatov, Sereda, 2001). It implies short period of time with the locally abnormal stress environment and rest of the time (up to present) with the normal stress environment around the existing displacement.
Some statistical models of hydrodynamic simulation of the present day fault-fracture system are based on representing of tortuous fluid flow patterns across sub-seismic fault population in average (Helth et al., 1994). Such stochastic models incorporate the mean length and mean density of the faults within the local calibration area within the basin (global scale). They could be treated as intermediate ones between the end-member fault transmissibility and fracture permeability models.
3.5. Chemical rock-fluid system response
The mineral changes at the molecular and crystal levels (diagenetic processes) as well as pore fluid changes on the solute component and phase content levels (pore water salinity and HC generation) represent mixed solid and fluid material response to basin processes indicated in the scheme in Fig. 5. The reason for incorporation of these essentially heterogeneous processes in a common group is that they all can be described within the range of the locally controlled mass conservation equation (8), where numerous parallel chemical reactions are integrated in current equilibrium of phases and solute components. The underground solutes are transported through porous and fractured media by Darcy flow, convection, dynamic dispersion and molecular diffusion (Verweij, 1993). Depending on driving forces the relevant transport equations at steady state assumption could be represented by generalized Darcy law (Darcy flow, convection) or Fick law (diffusion) or some mixed model like Sweeney's model (Sweeney et al., 1987) or Arrhenius Kinetic model (Tissot, Welte, 1978). In our target context, the chemical geo-fluid system response normally leads to more spatially local pore pressure variations in comparison to the mechanical rock response and PVT fluid response. Evidently, this higher frequency (spatially) process in both basin spatial and time scales seems to be less certain deterministically and consequently its commutative impact on present day pore pressure regime is less predictable. Still, some of them, such as oil and, especially, gas generation from source rock or local cementation of reservoirs are important origins of pressure generation and sealing mechanisms, which result in the high and sharp abnormal overpressure zones. Thus, it is important to recognise them among control parameters of the effective basin model aimed to be used for calibration and prediction purposes.
In our context, it is convenient to distinguish chemical between the responses occurring in highly and poorly permeable rocks. The first group of processes associated with the term "Cementation" controls the secondary porosity changes in terrigenous rocks and carbonates caused by solution O precipitation of solid material in mobile or immobile pore fluid. The cementation process is closely related to the pressure retaining mechanisms which reveal themselves via lateral and vertical pressure barriers formed by the sharp local reductions in rock permeability (tectonic and lithological seals) (Hunt, 1990; Weedman et al., 1998). The second group of processes associated with the term "Shale diagenesis" incorporates clay mineral and phase transform reactions described by Kinetic multi-reactions models with significant feed back. These additional fluid volume generation processes lead to a severe temporary excess of balance between available and required pore space in the representative rock volume (Barker, 1988; Dutta, 1985; Spencer, 1987; Sweeney et al., 1987).
The cementation processes in general are described by locally controlled conservation of mass equation (8) with the diffusion and conduction coefficients coming from relevant transport equations (10-12). Being considered more closely (Sweeney et al., 1987; Waples, Kamata, 1993), they can be subdivided on Pressure solution, Welding (Primary Cementation) and Secondary Cementation physical-chemical processes.
The pressure solution process is more important for carbonates than for sandstones. This term refers to the dissolution of mineral matter and the ensuing collapse of the mineral structure at a certain stress and temperature environment (Spiers et al., 1990), i.e. to effective porosity reduction. This diffusion-controlled process has a triggerlike mechanism with certain minimum depth or porosity levels to be reached by any given lithology before pressure solution starts (600-1000 m for limestone; 300 m for chalk; about 30 % - porosity for sandstone (Moore, 1989)).
The rate of pressure solution is described in the Fick law approximation (12) where the concentration gradient of diffusing ions is strongly dependant on specific grain-to-grain contact micro-area, which is in turn
a function of current stress and temperature. The simplified empirical relationships (Waples, Kamata, 1993) determine the specific volume of mineral matter mobilized during pressure solution by using the recursive multilevel model based on seven independent constants and two switches required to be specified/calibrated separately for limestone, chalk and sandstone.
The welding process is inverse to pressure solution. It implies precipitation of the solid material dissolved during pressure solution. It also can be treated as a porosity reduction due to welding of grains with the former type of cement (primary cementation). The pressure solution-welding can be considered as a locally conservative system where the amount of former cement material used for welding during single time step is limited by the mineral matter mass dissolved during Pressure Solution (Rittenhouse, 1971) within the open pore space. The relevant quantitative analysis of porosity reduction due to welding is based on the Porosity Solution model with an extension to account for linear welding rate increase in continuously reduced pore space. This extension requires one more independent lithology constant to be specified/calibrated.
The secondary cementation process is relatively more global than welding, since it implies transportation of dissolved material through porous media over a greater distance driven mainly by Darcy flow (10). The rate of cementation consequently is an inverse function of the fluid conductivity (11) of the host aquifer formation defined with dissolution-precipitation constants to be specified/calibrated for mineral matter of cement. These constants can be made available with acceptable accuracy based on petrographic observation and analysis generalized for the given basin environment as a whole. All together they allow definition of the diagenetic sequence of cementation-dissolution and provide chemical constraints on pressure seal built up. As an example, for the Gulf of Mexico the following diagenetic cementation line was reported (Weedman et al, 1998):
• T < 90°C - quartz-overgrowth precipitation
• T ~ 100°C - calcite cement precipitation
• T ~ 125°C - dolomite replacement of calcite; carbonate cement dissolution; chlorite rim cement and prismatic quartz precipitation
• T ~ 135°C - precipitation of pore filling kaolinite and calcite.
The effective porosity losses due to cementation are defined as an inverse to the hydraulic conductivity function extended for the Pressure dissolution rate and solubilities of calcite and silica fractions (Waples, Kamata, 1993). Last two lithology controlling constants are also a matter of specification/calibration.
Thus, all together pressure solution, welding and cementation processes can be considered as underlying the porosity reduction mechanisms in connection with chemical response of permeable fluid system. However, using such a set of parallel reactions with non-linear feedback links could be feasible for forward simulation but not for calibration of the relevant model.
A cumulative robust approach to modelling of secondary porosity and permeability reduction was developed by Schmoker and Hester. They suggested to use power function of the thermal maturity (M) to define the secondary porosity reducing trend against temperature scale (Rittenhouse, 1971) as the following:
where A and B - the empirical constants and M - the thermal maturity indicator represented by either TTI or vitrinite reflectance which depends exponentially on the temperature T and linearly on the basin time t (Rittenhouse, 1971).
By analogy with the flexure mechanical response controlling primary porosity reduction via effective consolidation constant C in global scale mass conservation equation (7*), the secondary porosity losses can be apparently controlled by effective bulk expansibility of the solid material with the temperature B. Again, depending on the level of details this constant could be treated as a cumulative EBM model parameter of the given formation or substituted explicitly by relevant combined empirical function (7*) or derived in average from the set of undergoing models introduced for parallel processes.
The terms "Clay" or "Shale diagenesis" are used to recognize the collection of many temperature-dependent reactions resulting in mineral and phase transforms. There is a long sequence of kinetically controlled first-order chemical processes. However, from the pressure generation viewpoint two of them are considered to be the most important: Smectite-Illite transformation (Dutta, 1985) and organic reach shale (Kerogen) degradation with the hydrocarbon generation (Sweeney, Burnham, 1990). As a series of several parallel first-order chemical reactions both processes are formally described by the same first-order differential model written for the basin time scale:
where m(t) - the current mass transform potential; k(t) - the constant of Arrenius defining the asymptotic rate of the reaction; E - the energy of activation of the reaction [kkal/Mol]; R - the gaseous constant; T(t) - the current temperature.
0= A [M (T, t)]B,
(19)
dm(t)/dt = - k(t)m(t), k(t) = A exp[-E/RT(t)],
(20) (20*)
The formal difference between clay mineral transformation and HC generation reduces to different activation energy6 and the relevant Arrenius constants. Physically, however, the first process reduces just to crystal level rearrangement with the releasing part of internal bound water and changing of pore space fabric, whereas the second process leads to solid - liquid - gaseous phase transforms with generation of HC fluid portions, which volumetrically exceed released pore space. Thus, the excess rate increases non-linearly when the gas generation temperature window is reached (Sweeney, Burnham, 1990). In connection with overpressure generation the effect of kinetically controlled clay diagenesis can be evaluated quantitatively (Gandy, Berg, 1997; Madatov, Sereda, 2000a). The general conclusion based on simulation of this process and data analysis (Swarbrick, Osborne, 1996) can be formulated as the following:
The contribution of Smectite-Illite transformation is estimated to be rather uncertain and modest among overpressuring mechanisms in sedimentary basins because the volume of fluid released is small and the dehydration is inhibited by the build-up of pressure.
The HC generation becomes primary significant overpressure mechanism when it reaches gas generation phase at about 120-140°C. Here, the specific volume requirements can exceed available in orders (Barker, 1990). The major problem with the gas generation/oil cracking is that the pressure buildup most likely retards the reaction (Osborne, Swarbrick, 1996).
The clay diagenesis response can be effectively accounted by introducing the relevant dynamic model for the source terms SQj and/or Ssj in both globally and/or locally controlled mass conservation equations (7-8), respectively.
3.6. PVT fluid response
This thermal and pressure controlled global scale process is entirely determined by the equations of state for pore fluid (13) in agreement with its fractional content. In contrast to the effective compaction and permeability trend functions defined empirically for the heterogeneous porous media in its effective solid material approximation (13* and 11) the equations of state for every given fluid component determine the dynamics of really homogeneous media. Thus, the PVT behaviour of the relevant fluid constants cj and fSj in globally controlled conservation equation (7*) can be predefined accurately for the given pore brine content.
In the general BM case the pore fluid should be represented by mineralised water occasionally with the immiscible and saturated hydrocarbon components (j = 2,J) in agreement with the current sj from (8) and its HC fraction content defined according to the source rock type and current phase balance (9). The solubility of HC in pore water is controlled by the saturation threshold diagram depending on the current temperature and pressure levels (Denesh, 1998; McCain, 1990). The physical properties of the pore fluid, which affect its density, are defined according to the "ideal mixture" law. In other words, it is the sum of the properties of the pure components, in relation to their mass fractions. The density of any "pure" phase is defined depending on the pressure and temperature via corresponding isobaric thermal expansibility ($) and isothermal compressibility (c) constants (see formula 7*).
In the pore pressure prediction context it is important to be able to estimate the density difference between generated hydrocarbons mixture and source rock (effective pressure source due to the primary HC migration) or between immiscible hydrocarbon cap and host pore brine (buoyancy pressure effect due to the secondary HC migration).
For the typical range of pressure analysis in the sedimentary basin (depth: 0-6 km, temperature: 10-200°C, pore pressure 0.1-50 MPa) the following linear fluid density model is accepted to be workable
(Domenico, Palciauscas, 1979):
Pj (P,T) = PJ (P0,T0)[1 + ¡3J (T - T0) + cj (P - P0)], (21)
where P0 and T0 - some reference "normal" conditions for pore pressure and temperature, where the control lab measurements are available; ¡3-the isobaric thermal expansibility, c - the isothermal compressibility.
For the "pure" gaseous phase of pore fluid the state equation of the ideal gas can be used to define its density versus P and T. In practice, some correction for "non ideal" gas behaviour under pressure must be made - (Z - factor, or Kats constant - X in formula (22).
P2 (P,T) = P2 (P0T0) [(P/P0)(Tc/A7)]. (22)
See also (Oil and Gas., 1984) for more details.
In nature, however, "pure" phases of pore fluid are not observed. For example, pore water becomes more and more mineralized with depth, so its density should gradually increase. This effect is partially compensated for by temperature extension and, as a result, the linear model (21) appears to be even more precise than a parabola (Madatov, Sereda, 2000a).
6 E is about 20 kcal/mole for Smectite-Illite transformation and about 55 kcal/mole for the Kerogen degradation.
Thus, the current PVT response of pore fluid in terms of its density and viscosity can be component-by-component controlled by the equation of state (13) with the isobaric extension and isothermal compressibility constants predefined in advance (McCain, 1990). The HC fractional content at the primary migration (generation) phase is less certain. Some Kerogen classifications (Tissot, Welte, 1978) or calibration based on inverting of the generation model (20-20*) (Fox, 1995) can be applied here. Finally, prediction of HC fractional content and cap elevation of the relevant immiscible fluid part at any given location requires full scale inversion of the relevant dynamic indicators with respect to the 3-D Basin Model. Based on the present day state of the BM approach it is hard to imagine a solution of this problem in the nearest future. Thus, the evaluation of buoyancy pressure effect could be only based on a priori given HC fractional content and cap elevation.
4. The effective HD basin model against the general Basin Model: Assumptions, restrictions, validity
The full-scale forward basin model described above appears to be redundant in the context of the model acceptable for 3-D simulation of overpressure phenomenon and further implementation for its prediction.
It is clear intuitively that the present day overpressure phenomenon is non-sensitive or weak sensitive to the most of local level mass conservation effects. For example, seven porosity reduction mechanisms considered in the context on mechanical and chemical rock response on stress and temperature increase during burial history (see s/sect.3.4-5) have twenty seven independent constants to be calibrated (Waples, Kamata, 1993). The amount of the porosity measurements and/or accuracy of the well-log to porosity transformations generally is not sufficient to support uniqueness and stability of the relevant data fitting problem. Indeed, such local level effects act simultaneously in nature and create cumulative effect of porosity losses. These effects cannot be uniquely recognised separately in retrospection of basin time scale based on the present day porosity data, i.e. without paleo markers of the relevant diagenetic processes. On the other hand, the available inverse overpressure models base on compaction trends and do not require distinguishing between these local level effects. Thus, the trend model approximated porosity losses in agreement with the given lithology seemed to be effectively applicable for the purposes. The effective approach to multi-mechanism porosity losses' process is illustrated in Fig. 8.
The review of the empirical background available for setting and calibration of the aimed effective basin model allows us to introduce certain assumptions in order to simplify the forward problem and define the class of the appropriate effective models. Inside of a certain class, the only real data amount and quality can serve as an objective constrain in further Earth model simplification in terms of upscaling. The formal way of such simplification consists in upscaling of the raw Earth model for overpressure simulation-prediction purposes. It has been considered in (Madatov, Sereda, 2003). Here we are going to consider the set of upper level assumptions
0,0
0,1
Porosity
0,2 0,3
Porosity losses mechanisms:
8,4
0,5
1000 ■
2000 ■
3000 ■
QL
a> Q
4000 ■
5000 ■
eooo
] 1 3 | □ —
□ □ J? f U
" «Gems nation 1 / J f a 1
□ □ □ □ □ □ J^r D r □ □ □ □
□ □ □ □ □ □
□ D :Mec nanical ;ompac
- □ □
□
□ □ □
TO
o c
TO -C
o u
TO
0
1 u
-C
o
SH SS CB
Repacking + + +
Crushing - + +
Ductile flow - + +
Pressure solution - - +
Welding - + +
Cementation - + +
Clay diagenesis + - -
0,0 0,1 0,2; 0,3 0,4 0,5
—=— Stress depended losses —■—Temperature depended losses
-Integra! losses - combine curve
-Depth trend approximation of combine curve
Fig. 8. Effective exponential approximation of combine mechanical and chemical mechanisms of porosity losses suitable for EBM of mixture lithology unit (speculative numerical example).
SH - shale; SS - sandstone;
CB - carbonate
which allows introduction of the class of the effectively invertable forward hydro-dynamic basin scale models (HD models) suitable for 3-D simulation of overpressure phenomenon and further implementation for its prediction.
Strictly speaking, every assumption requires formal and/or experimental (numerical modelling) grounds to be accepted. However, this way of argumentation seems to be far beyond the framework of the single paper. Instead, we will continue to follow the principles of maximal pragmatism and effectiveness of the framework forward model introduced in the first year report (Madatov, Sereda, 2003). Namely:
A properly upscaled purposely built Earth model must reproduce predefined target data set by its synthetic prototype with a misfit quality not less than fixed in advance error level. The optimally simple model parameter space is one, which has minimum orthogonal model parameters to be calibrated.
The reviewed above inverse models can be considered as implementations of the pragmatic approach to the pore pressure prediction purposes at rather rough static level of the basin model concept. Failures of this approach in more complex geological conditions than those where it has been originally established are due to ignoring the fluid dynamics aspects in simplified static models. The constant link with the general BM equations and principles of maximal pragmatism are prerequisite of successive accommodation of the forward basin scale models to more extended requirements coming from practice of pore pressure prediction.
Let us now consider some basic assumptions which must be done in order to reduce general basin model equations (6-9) to the form of combine 3-D HD model suitable for the fast and real data consistent overpressure modelling.
First, the set of physical assumptions is introduced. Then two more assumptions are discussing in connection with the numerical solution (see s/section 4.3).
4.1. Physical assumptions
1. Density assumptions
These assumptions imply that the changes in bulk density of rock are only due to porosity losses, which include both mechanical (compaction) and chemical (diagenesis) effects.
Thus, the changes in solid single mineral phase density versus pressure and temperature are ignored in mass conservation balance, unless some phase transform occurs. In the last case, the sample volume U (see Fig. 3) is getting pore fluid source and solid phase sink terms equal by mass. In other words, the terms SQ0 and ¿j2i in (7) are equal but have an inverse sign. Practically, this assumption allows introduction and specification of a source term as a lithology depended model parameter responsible for possible during loading phase transformations (Madatov, Sereda, 2000a).
Here, the term "Source" is not obligatory associated with the organic reach source rocks (Kerogen). The inside-crystal water release during clay mineral diagenesis (Smectite ^ Illite transformation, for example) can also represent an additional fluid generation taking place at certain stress-temperature conditions (Osborne, Swabrick, 1997). Thus, both mechanisms are potentially volume generating. Moreover, as it was shown in s/sect. 3.5, similar rock mineral transformations can simulate based on the same first-order kinetic models (Dutta, 1985; Sweeney, Burnham, 1990).
Thus, the kinetically controlled mechanisms of rock mineral conversion are effectively accounted in
global level mass balance equation (7*) by the effective fluid generation term SG1 defined as the following:
t
3G1(t) = A(1 - 0(A>/A - 1)'k(t)exp[- | k(t) dt], (23)
0
where A - the mass source rock potential, ^ - the porosity; p1, p0 - the solid phase and generated geo-fluid densities, respectively; k(t) - the Arrenius constant defining the asymptotic rate of the reaction (at infinitely high temperatures) in [1/time] units.
The derivation of this formula and numerical simulation of the relevant source term are available in (Madatov, Sereda, 2000a).
The changes in pore fluid densities are expressible in linearised form (7**) (Domenico, Palciauscas, 1979). The contribution of each phase is controlled by the PVT equation of state (Denesh, 1998). See also Appendix 2 and (Madatov, Sereda, 2000a) for more details.
2. Stress assumptions
The effective stress assumed to be defined from Terzaghi's model (2) or (5). Hence, the solid phase movement relative to global (stationer) coordinate system has only the vertical non-zero component and upper equation in (7) can be reduced to the form:
d[f0(1 - 0]/dt = A(1 - 0 q0 - Q (24)
It is important to stress, that the empirical relationships for effective rheology and underground temperature variability are practically available in the form of Normal Compaction depth trends for each given lithology type and thermal depth gradient for the given basin area, respectively.
3. Compaction assumptions
The changes in the specific volume (porosity) is accounted by Athy formula - <j) = exp(-rZ) (Magara, 1978) that incorporates the primary (compaction) and secondary (diagenesis) porosity losses as the following:
r= C [Vz(c)] + B [Vz(T)], (25)
where C, B - the bulk compressibility and expansibility of the porous media; Vz(o), Vz(T) - the depth gradient of the stress and temperature increases, respectively.
Thus, the four terms' model (7***) is substituted by twice simpler expression:
(1/^)(d^/dt) = -r(dZ /dt), (26)
where the instant burial speed could be effectively treated as an current sedimentation rate.
4. Single phase assumptions
The multi-phase pore fluid assumed to be effectively replaced with the single-phase one (composite liquid). The effectiveness of such replacement is defined by account in (7) for changes in composite density according to the thermodynamic law (PVT properties) but ignoring multi-phase fluid dynamics given by (8). Thus, the term "fluid dynamics" could be effectively treated as a hydrodynamics (HD) in this single-phase approximation.
5. Heat transfer assumptions
The temperature variation can be approximated by the thermal gradient for the given time-space interval of burial history. Thus, instead of solving heat conservation problem (9) at each time step of burial history the common for the area temperature gradient is used in (7*) to represent the temperature derivative as the following:
(dT/dt) = grad, (T)(dZ/dt)j, (27)
where the instant burial speed at the j-th episode of burial history could be effectively treated as an current sedimentation rate.
The convection term in Darcy law (10) is effectively accounted by PVT controlled fluid density variation (see assumption 1). Thus, the relative speed of steady state fluid flow is approximating by the linear function of the pressure gradient only:
q = DVP, (28)
where conduction coefficient D is defined in single phase (hydro-dynamic) fluid and isotropic rock matrix approximation of general form (11). In other words, K0 assumed to be a scalar coefficient and Kj = 1 in formula (11).
6. Effective 3-D fluid flow assumptions
The fact of different pressure regimes observable in well permeable (aquifer) and poor permeable (aquitard) intervals of the section associated with the sand and clay lithology, respectively, is well proven (England et al., 1987; Stump et al., 1998; Bredehoeft et al., 1998). The key factor underlying this difference is rather different rate of the excess pressure equalisation and direction of the relevant fluid flow. In particular, the lateral component of fluid flow during releasing pore water from compacting host rock is significant for sand intervals whereas it is negligible for clay-shale intervals. Thus, there is no need to resolve the 3-D HD basin scale problem for whole fluid system but rather for the aquifer part of it, which effectively works as a pressure communicator during burial history.
The keystone here is assumption that 3-D Earth model during its burial history can keep this feature preserved. The HD system relevant to this assumption allows to approximate description in terms of finite set of aquifer formations sandwiched within the host aquitard formations. The fluid flow divergence term in (7*) for such HD system is representable in the form of additive explicit and implicit sink components according to the following considerations.
Firstly, the divergence term is subdivided on the explicit vertical sink term - d(DZ (cP/dz))/dz and the implicit lateral sink term - sqL as the following:
div(q) = d/dz [Dz (dP/dz)] + SqL. (29)
The relevant 1.5-D HD model (Madatov et al., 1996) describes discharge of instantly charged 1-D across formation model relevant to the 1-D burial history concept (see s/sect. 3.2). All 3-D pressure communication effects in this model are effectively represented by the implicit lateral sink term - 8qL. Secondly, the same term is represented dissymetrically to (29) in the form:
div(q) = d/dx [Dx(dP/dx)] + d/dy [Dy (dP/dy)] + SqV, (29*)
where lateral HD conductivity components Dx and DY are fully controlled by the aquifer formation properties; 8qV - the implicit vertical sink term.
The relevant 2.5-D HD model (Madatov, Sereda, 2001) describes discharge of instantly charged 2-D aquifer formation model built based on 2-D interpolation between 1.5-D HD models described by (29).
The relative contribution of lateral sink term in mass balance for the aquitard formation assumed to be negligibly small in comparison to the amount of fluid releasing in the direction of z axis (England et al, 1987; Bredehoeft et al., 1998; Madatov, Sereda, 2000b). The rest of the section, represented by the aquifer formations, provides an effective lateral pressure communication of the relevant 3-D HD system during burial history. Consequent use of equations (29) and (29*) allows effective 3-D HD forward modelling for the relevant Earth model.
The more detailed consideration of this combine model with some computer examples attached are given in the next sub-section.
4.2. Basis equation in reduced form
The implementation of the listed above basic assumptions allows reducing of basis globally controlled mass conservation equations (7-7*) to the following upscaled 1.5-D form:
Aj (dP/dt) = {Bj [grad(7)]j + CJ}(dZ/dt)J + 5q + div (q), { div(q) = d/dx[DXj (dP/dx)] + d/dy[DYj (dP/dy)] + 8qV] Vj e {Uaf}, (30)
div(q) = SqLJ + d/dZ [DZj (dP/dZ)] Vj g {Uaf},
where the sub-script index j instead of "0" for the solid phase used in the above formulas (7-7*; 23-24) indicates now the formation unit number within the consequent set of sedimentation history events {U} and the sub-set {Uaf}c{ U} indicates a sequence of the aquifer formation units recognised among others. Due to single phase assumption 4 (HD model) we omit now the sub-script index "1" at geo-fluid state constants.
Included in (30) terms essentially represent two groups of the approximate descriptors: functions of the geo-fluid system response - {Aj, B, Cj, Dj, Gj }; functions of the framework basin process - {(dZ/dt) [grad(T)]j}.
They both can be interpreted in the agreement with background empirical models introduced in s/sect. 2.1. In particular:
Aj is the combine matrix-fluid compressibility term in [1/Pa] given by:
Aj = [Cj% -p)g - 4c\. (31)
Cj is the matrix compaction term in [1/m] given by (see assumption 3), c - the isothermal porous water compressibility constant in [1/Pa]:
Cj = j(1 - 0. (32)
B - the geo-fluid thermal expansibility term in [1/C] independent on the formation episode number and given by:
B = <I>P, (33)
where ¡3- the isobaric porous water expansibility constant in [1/C].
Dj (DXj, DYj,DZj) - the hydraulic conductivity vector in [D/(sPa)]7 given by (see assumption 4 and formula (11):
Dj(z) = Kj [S;^(z)Mz), (34)
where Kj - the lithology dependent rock matrix permeability vector; n - the viscosity of pore water vs. depth function. The rock matrix permeability term is getting from established for the specific lithology trend function of current porosity (see s/sect. 3.4 and formulas (17-17***) for more details). The lateral components of aquifer formation permeability affected by faults are mainly controlled by fault permeability (Helth et al., 1994; Antonellini, Aydin, 1995; Manzocchi et al., 1998; Madatov, Sereda, 2001). The pore water viscosity is defined as a function of current depth at the given temperature gradient in agreement with available empirical relationships (for example, according to the Mercer's formula (Denesh, 1998)) and assumption 5. SGij - the effective dimensionless fluid generation term given by formula (23) (see assumption 1).
The meaning of terms {(dZ/dt), [grad(7)]j} directly follows from assumptions 5 and 6, respectively. These processes are also discussing in s/sect. 3.2 and 3.3.
In addition, there are some assumptions, which implicitly follow from the final difference method of the numerical solution of the relevant forward problem (30).
4.3. Quantization assumptions
Some important assumptions must be mentioned in the context of the numerical solutions of the forward problem (30). The general scheme and numerical features of the relevant combine 3-D forward modelling operator is discussing in the next sub-section.
7 It is more suitable for basin scale HD modeling to use the following conversion into basin time scale expressed in million years before present time (mY): 1 [D/(Paxs)] = 0.317 [m2/(PaxmY)].
7. The 1-D burial history assumptions
The stress & temperature evolutions in agreement with assumptions 2 and 4 imply discrete dynamic models with steady load & heat increment during the finite number of sedimentation history episodes (the number of formation units).
If so, the relevant burial history for any single 1.5-D components of the combine 3-D HD basin model could be uniquely recovered by using backstripping technique (Allen, Allen, 1990). The burial history for such model is subdividing on the relevant number of sedimentation events - macro time steps associated with the constant sedimentation, heating and lateral drainage rates.
Evidently, the positive burial history events (macro steps) are associated with sediment depositing - positive increment in loading (assumption 2), heating (assumption 5) and further compaction of sediments (assumption 3). Conversely, the negative burial history events are associated with erosion - negative increment in loading (assumption 2), cooling (assumption 5) and keeping porosity at minimal reached before level (assumption 3).
8. The charge-discharge cycle assumptions
The positive and negative load events (macro steps) in the discrete 1-D approximation of the burial history can be associated with the positive (charge) and negative (discharge) cycles in discrete approximation of the relevant HD system. The finite fluid potential collected during every j-th charge cycle can be defined based on system response functions {Aj, Cj, Dj, Gj} in agreement with reached burial and temperature rates according to framework basin process assumptions 3 and 5.
The discharge of HD basin system occurs during the same time cycle in agreement with combine 3-D steady state fluid release model (29-29*) in agreement with assumption 6.
The single charge-discharge macro cycle can be described by final difference approximation of the relevant HD forward problem (Transient Flow problem) (Wang, Anderson, 1982).
Listed above assumptions reveal the main advantage of the effective basin model (EBM) and relevant numerical modelling against the general BM concept, which is possibility to use practically for data inverting. In agreement with the overpressure phenomena context, this is the most important prerequisite of successful solution of the prediction problems.
The unavoidable price for this advantage is the restriction of the relevant EBM validity.
4.4. Restrictions
The general restrictions of the introduced effective basin model follows from the above described assumptions (1-8). In the general context of the reviewed BM approach, they can reduce to the following items.
1. The applicability of the EBM concept for pore pressure prediction is restricted with the Extensional type of sedimentary basins (see s/sect. 3.1). The influence of any lateral compression components on local stress environment (non-vertical fault or intrusion activation) is not included explicitly into the EBM in terms of relevant independent model parameters.
2. The thermal history is approximated with the steady heat flow through homogeneous media. Such simplification ignores any local spatial variation of the basin temperature regime associated with local heat sources movements, abnormal thermal conductivity zones, etc. (see s/sect. 3.3).
3. The mechanical (primary) and chemical (secondary) compaction mechanisms are not distinguished in Normal Compaction Trend approximation for the given lithology. In other words, any local level environmental porosity variations are excluded from predefined. Thus, relevant forward modelling of overpressure evolution based on the given EBM appears to be not sensitive to the local level primary and secondary porosity reduction mechanisms.
4. The 1-D burial history does not account for 3-D sedimentation phenomena like non-compensated deposition, turbid flow deposition, erosion of dipping formation beds, fault and salt tectonics, etc. As a result, composite 3-D formation models cannot account for apparent disappearance / appearance of anticipated formation units due to intersecting of the model with a deviated fault path, apparent disappearance / appearance of salt layer associated with distortion of anticipated formation thickness.
5. The effects of HC (gaseous, oil) phase generation and secondary migrations are accounted in a single phase approximation. As a result the relevant control parameter is not consistent with generally accepted for source rock characteristics like TOC or HC Index. The buoyancy and/or vapour lock effects (Swabrick et al, 1999) on overpressuring can only be accounted effectively.
These restrictions shrink the applicability of the discussed EBM concept for overpressure prediction to certain global level geology setting and reasonable assumptions about acting overpressure mechanisms. Still, in comparison with reviewed above inverse models widespread in the exploration engineering the EBM class provides much more general and universal ground for Earth model setting.
In conclusion of this sub-section it is important to pick out, that consequent constraining of the general basin model with purposely oriented assumptions and substitution of multi-level local scale problems with the effective single-
level inverse model allows approaching to the practically invertable level of the aimed overpressure model. In agreement with the introduced above principle of Maximal Pragmatism all diversity of geological processes and relevant geo-fluid system responses leading to well proven overpressuring mechanisms (see Fig. 5) can now be effectively substituted with the EBM given by the combine 3-D HD model. The relevant forward modelling operator appears to be determined on a few-dimensional model parameter space, which includes three formation parameters controlling compaction (31-32) HD conduction (34) and volume generation processes (23) plus one fault transmissibility parameter controlling lateral pressure communication in aquifer formation (18*). These parameters must be specified in terms of most likely values and variation range for each given lithology types and fault morphology (Madatov, Sereda, 2000a; 2001). The relevant set of control parameters is a matter of sensitivity analysis and model calibration (Madatov, Sereda, 2001).
5. The 3-D combine overpressure solution
Strictly speaking, neither sub-vertically (29) nor laterally (29*) charging-discharging HD models of pore pressure evolution provide physically grounded interpretation of relevant sink and recharge terms. However, they jointly can supplement each other up to establishing an effective and physically valid 3-D HD model in the basin scale. The most attractive feature of this combined model is its simplicity in terms of model space dimension (Madatov, Sereda, 2000a), required cells number to be specified (Madatov, Sereda, 2001) and computing speed. As it was stated above, this is a crucially important prerequisite of the HD model applicability for calibration and for pore pressure prediction purposes.
There are many proofs of the fluid transport redistribution between aquifer formations and bounding shale rocks which occurs during sedimentation (England et al., 1987; Stump et al, 1998). The best examples of rather different excess pressure regimes inside of aquitars and aquifers are related with well known Centroid effects (Traugott, 1996) and pressure compartmentalisation (Leonard, 1993; Yardley, 1998).
It is important to stress, that not all the depth intervals of the present day geological section but just their small part provides a real lateral HD connectivity during burial history. This is especially important for the oldest part of the sections overlapped with the thick poor permeable mud-rocks and carbonates. As it was reported from the North Sea case study (England et al., 1987), the only small fraction of whole section in thickness representation (about 10 %) provides main part (about 90 %) of the total volumetric lateral outflow from the whole sediment column during its compaction. The similar conclusion came from the Caspian Sea case study (Bredehoeft et al., 1998). Our own modelling and calibration results achieved in the North Sea, Pechora Sea, Caspian Sea and Western Canadian offshore zones also support this conclusion (Madatov, Sereda, 2000b).
The forward modelling results published in 2001 (Madatov, Sereda, 2001) also reveal that the volumetric ratio between the portions of pore fluid releasing in the vertical and lateral directions differs for aquifer and aquitard formations in orders. Namely, the sink terms power in steady state approximation could be estimated as 10-2-10 2 sp. volume units per 1 mY for good permeable sandy aquifers and less than 10-9sp. volume units per 1 mY for bounding shale rocks.
In the context of the aimed computer simulation of HD basin system, this fact means nothing else, but possibility to treat a bigger part of the sedimentary rock section as an essentially 1-D charging-discharging model introduced above. Namely, a poor permeable and fine grained sediments (aquitards) provide mostly subvertical directions for pore water release during compaction. Naturally, high permeable sediments (aquifers) are mainly controlling the lateral connectivity of the whole section. Thus, 3-D charging-discharging HD modelling can mainly base on aquifer geometry and conductivity evolution through the burial history time, where enclosing aquitards supply them with the additional portions of the pore water to be released laterally.
Below the important aspects of numerical solution described in more details and some results are discussed.
5.1. General scheme
We will consider now the general scheme suggested for the effective forward modelling and calibration of the combined HD basin system. It will be based on consequent application of vertically and laterally discharging models of fluid flow divergence (29) and (29*) to the same earth model.
Let the set of individual 1.5-D HD models be established within the area. Let also each j-th of the individual 1.5-D HD model include the same aquifer formation as a common aquifer unit providing lateral pressure communication for the relevant HD geo-fluid system (see Fig. 9). Clearly, the vertical component of the fluid flow coming through the top and bottom of this aquifer unit at any offset well location can be reconstructed by any time slice of burial history based on calibrated geo-fluid system response - {Aj,Cj, Dj, Gj} of the relevant j-th 1.5-D HD model8 (see descriptors to formula 30). Hence, the vertical part of the fluid divergence (the vertical recharge term SZqAU) can be easily calculated. At that, the assumption about absence of additional fluid generation inside of the aquifer formation must be done during every time step. An interpolation technique can be implemented to estimate
8 In terms of the model calibration this set of responses at predefined empirically functional links with the relevant coefficients represents model parameter vectors {xb x^ x3, x4}j sorted out according to formation unit consequence.
a
b
Fig. 9. Incorporating the set of 1.5-D single HD models calibrated against offset wells (Framework Pressure Model) with 2-D grid model of aquifer formation in order to get a combined 3-D solution. a - the 3-D view of the aquifer formation penetrated with the set of calibration well (red lines);
b - the structure map of aquifer unit (AU) with well
mark on it (red circles); c - the cross section along green trace with the combine sample and the relevant fluid flow patterns highlighted
the AU parameters including the vertical sink term SZqAU required for 2.5-D modelling at any intermediate x, y position between the offset wells. It can be done by any slice of the basin history time grid.
The forward problem (30) for the given aquifer formation now can be numerically resolved on the relevant spatial grid. By analogy with the vertical component of the fluid divergence, its lateral part SLqAU can also be easily calculated for the given aquifer formation at any of time steps. In contrast to the effective meaning of the lateral sink term at the 1.5-D HD modelling, its estimation now will be based on the explicit solution of 2.5-D lateral problem (29*).
Now the updated and corrected for 2-D HD communication lateral sink terms can be returned back and redistributed among the set of individual 1.5-D models according to their x, y positions on corresponding 2-D grid (see Fig. 9a and c). After this stage the calculation could be repeated from the start point.
The flow diagram of the relevant combine algorithm is presented in Fig. 10.
Fig. 10. The mutual exchanging of the 1.5-D HD across formation (above, blue) and 2.5-D lateral (below, green) forward problem solutions at the combine 3-D HD modelling.
^.5-0/2.5-d - the forward HD modelling operator applicable for 1-D and 2-D Earth models, respectively; SqZ,ScjL - the vertical and lateral sink terms, respectively; f - the adjustable formation and fault parameters, respectively.
5.2. Numerical solution
The numerical solution of the forward problem (30) in its final difference approximation for 1.5-D case (29) was considered in details in (Madatov, Sereda, 2000a). The absolutely stable implicit scheme of 2.5-D case (29*) was introduced in (Madatov, Sereda, 2001).
Let us now consider combining of 1.5-D and 2.5-D solutions of the relevant 3-D HD forward problem more rigorously. To do that, one must formulate it together with boundary and initial conditions written for 3-D space in the basin time scale.
The general mass conservation problem (30) will be posed now within the 3-D area Q, for basin time interval [T0,TM]. Let the beginning of burial history (T0) be associated with the oldest formation unit described in area stratigraphy column and its end (TM) - with the present time. The forward solution of this HD problem with respect to excess pore pressure is searching inside the 3-D target area Q given by:
{
y™< y < y» (35)
0 < z < zmax.
Since the forward problem (30) was formulated above for the hydrodynamic system (assumption 4), the pore fluid considered through all time-space intervals of the basin model is assumed to be pore water with all its PVT properties (Denesh, 1998).
The terms Aj, Cj, Dj, Gj of basis equation (30) are essentially functions of lithology. According to the specification (31, 32, 34 and 23), they are uniquely linked with each given formation unit sedimentation history through the relevant model parameters: compaction constants (26), specific surface constant (17) and mass source rock potential (23). It means that the relevant components of model parameter vector {xb x2, x3, x4}j could uniquely be determined for every given sample volume Uj of loading/unloading formation rock at any time slice if its paleo (x, y, z) position (see Fig. 3) is recoverable and the relevant empirical relationships or PVT laws are predefined.
The uniqueness in restoring of sedimentation history in terms of paleo geometry, paleo stress, paleo temperature is a matter of data quality available for description of present day geology, geophysics measurements and calibration technique. Full scale investigation of this matter is far beyond the discussing problems. In connection with actual aims of the given discussion, one could stay with pragmatic logic described above. From this point of view the assumptions about terms Aj, Cj, Dj, Gj are acceptable for developing suitable approach to the basin HD system modelling-calibration-prediction.
In particular, the rock matrix compressibility (A), the rock matrix compaction (C) and HD conductivity (D) for most of the clastic and carbonate sedimentary rocks can be predicted in advance based on "Normal Trend" relationships (Magara, 1978). If so, the functions A(x,y,z,t), C(x,y,z,t) and D(x,y,z,t) are defined with accuracy restricted by validity of empirical relationships and quality of relevant model parameters' calibration. The cumulative source term function G (x,y,z,t) is generally more uncertain. Still in tectonically relaxed basins this function can also be defined via subsidence, paleo water depth, paleo temperature gradient curves (Allen, Allen, 1990).
Strictly speaking, the mass balance between release and generation of pore fluid for finite 3-D area (35) cannot be assured without suitable boundary conditions.
Let G be now a full set of surfaces staking the spatial boundary of our target 3-D area. We can determine G in the following form:
G = Gt U Go U G, U Gb, (36)
where:
Gt - is the upper surface (erosion basis) where excess pore pressure is stationary equal to zero according to assumption 1.
Gb - is the bottom part of the 3-D area coinciding with the current position of the non-permeable basement.
Thus, a non-flow part of the 3-D area can formally be associated with this surface. Go - is the part of side surfaces that provides connection of the 3-D area with the Remote Discharge Zone via available aquifer formations (Madatov, Sereda, 2001). In agreement with this the lateral excess pressure gradient for the relevant edge cells is determining via the current pressure value and offset distance value, which is constant for the area. Gi - is the part of side surfaces that provides connection of the 3-D area with the neighbouring areas via the steady state pore water flow regime.
The boundary conditions for the forward problem (30) with boundary architecture (35-36) now look as follows:
{P (x, y, z, t) = 0, V(x,y,z) e Gt,
P (x, y, z, t) = 0, V(x,y,z) e Gt, (37)
d(P (x, y, z, t)) /cn = P(x, y, z, t)/A1, V(x,y,z) e G0,
P (x, y, z, t) = P0(x,y,z), V(x,y,z) e Gi.
The initial conditions by the start of the burial history model assumed to be known in the form:
P(x,y,z,0) | (xy,z)e n = P0 (x,y,z,) | (xy,z)e fi. (38)
Practically the excess pore pressure by the start of the oldest formation depositing supposed to be equal to zero.
Note, that neither Gt nor Gb surfaces remain spatially fixed during the interval [T0,TM] of burial history. Their positions at any te [T0,TM] are determining according to the paleo water depth curve and subsidence curve (Guidish et al., 1985). In connection with the numerical solution characteristics, the problem (35) is posed for the 3-D area with the time depended z grid size for every x, y cell.
Finally, the combine algorithm of 3-D forward problem solution incorporating consequent 1.5-D and 2.5-D models can represent as the following.
Let the regular 3-D grid ^?x,y,z to be reduced down to optimal irregular grid coTxy^>diz according to upscaling strategy (Madatov, Sereda, 2003). This grid is composed from a certain number MV of 1.5-D HD models given on the
relevant 1-D grid oRz and ML number of 2.5-D HD aquifer formation models given on 2-D mesh coTxr At that ML aquifer units provide a lateral pressure communication within the target area through all its burial history according to their age -position on time interval [0,7] = [T0,TM]. Note, that MV nodes of triangular mesh are used for reproducing the AU geometry. The relevant 1-D across formation models tied to nodes are specified with the optimal number R of formation units established in the area stratigraphy column (see Fig. 8). The set of model parameters required for 1.5-D and 2.5-D HD modelling can be optimally determined based on the same upscaling strategy. At that, the model parameter set for 1.5-D and 2.5-D operators require specification of lateral and vertical sink terms, respectively (Fig. 10).
In agreement with the formation unit definition given in section 2 (see Fig. 3), the full time history interval can be subdivided on M stratigraphy intervals, where the sedimentation, subsidence and sea bottom variation rates in average assumed to be constant. The consequence of the relevant macro time steps covers the full time gate as follows: M
[0,T] = u [Tm-1, Tm], T0 = 0, Tm = T. (39)
m =1
The combined modelling approach implies implicit routine process, in which the solutions currently obtained for 1.5-D and 2.5-D forward problems on the regular time grid cJRt inside of any macro time step are linked on the common spatial grid coTx,y^>dRz and are currently exchanging by sink terms as it is shown in Fig. 10.
We will consider now an arbitrary m-th macro time step from the list (39) in order to get step by step description of this process. Let by the end of (m - 1)-th macro time step the excess pressure solution Px,y,zm-1 be available for any nodes of the spatial grid coTx,y^>dRz together with distribution of formation units and associated model parameters. This pressure solution can serve as an initial condition (38) by the start of next m-th macro time step. The regular local time grid ofit contains constant for every macro time steps number NS of the simple time steps ht and Nc of the single charging-discharging cycle hT. At that, the pair step-cycle must satisfy condition: hT / ht = l, where l > 1 is integer. The length of the single charging-discharging cycle is choosing to get an instant charge term at any x, y position of the 1-D model in oJxy grid according to the quantization rule (assumption 8). This average pseudo source term is used to simulate the cumulative excess pressure impact during a single tact of the instant charging of 1.5-D HD system in agreement with the basin process rate.
Let the values of the single time step and charging-discharging cycle be fixed. The combined solution of 3-D problem (30) now includes the following repeating steps:
STEP 1. Charge instantly each 1.5-D HD model of the Mv set according to the relevant averaged charge cycle. STEP 2. Get Mv individual solutions9 Px,y,z(1/2) on half of charging-discharging cycle distributed within the spatial mesh coTx,y.
STEP 3. Collect the current value of A(x,y,z), C(x,y,z), D(x,y,z) and SqZ(x,y,z)w for every 1-D models tied to ML
aquifer units and interpolate them on the relevant regular grids cJRx,y attached. STEP 4. Get ML individual solutions11 Px,y,z(1) on the second half of charging-discharging cycle distributed within
the set of regular grids {fiRRx,y}. STEP 5. Redistribute solutions Pxy,z(X> and ¿qL(x,y,z)12 back to the irregular grid coTxy^dRz according to x,y position of
the 1-D models tied to ML aquifer units. STEP 6. Make a single increment in the charge-discharge cycle. If the corresponding time interval Tm related with m-th formation unit depositing is not over - go to step 1. Else - start from step 1 the m+1-th charge-discharge loop for new formation unit.
The six-step routine is repeating M times. The excess pressure solution Px,yJM) by the end of the given macro time step is available on the optimal spatial grid coTxy^>dRz upscaled for calibration purposes. Thus, the set
of 3-D pressure models {P^1, P^2,.......PJ1, ... Pa1} effectively available for the target area Q, covers the full
time interval [0,T] of the burial history. Clearly, every intermediate solution represents distribution of paleo-excess pressure, whereas the final one reproduces the present day picture and can be used for calibrating against available formation pressure data.
Some important numerical features of the described combined solution must be highlighted.
Stability
The term "stability" is used here to indicate the reliability of the numerical solution of the forward HD problem wherever the model parameter vector x is defined inside of the physically grounded range of its possible
9 Every solution PzyT) obtained for the depth grid and attached to the given x, y 1-D model position is based on the ADI algorithm (Madatov, Sereda, 2000a).
10 The instant vertical sink term can be made available based on the current 1.5-D pressure solution (Madatov, Sereda, 2000a).
11 Every solution PxJl) obtained for the spatial grid aRxy and attached to the given AU position is based on the ADI algorithm described in Appendix 1.
12 The instant lateral sink term can be made available based on the current 2.5-D pressure solution for the given aquifer unit.
variation. Note that in general, the numerical implicit scheme based on the ADI method provides an absolute stability (Samarsky, Gulin, 1989). It means that any small perturbation of the model parameter vector results to a finite misfit between reference and testing synthetic values. At that, however, the properties of time and spatial grids can influence on the final result in terms of numerical modelling errors E(x). It is important to guaranty that the average level of E(x) value is invariant to the specified position of the model vector x within a priori given subset Xa of the model parameter space X. As soon as E(x) vs. grid parameters13 can be treated as a stationary process it becomes property of a chosen final difference scheme, not a model property. In this case, the combined 3-D HD modelling operator determined jointly with its grid property can be treated as a numerically stable one.
The results of the numerical test are summarised in Fig. 11. The test was built as follows.
First, the 3-D HD model x was created based on the pressure compartmentalisation example (section 6).
Second, the synthetic prototype M"J':i'(x) was computed based on the developed combine scheme with the operator defined on the knowingly dense grid.
Finally, the time grid step was gradually increased and the misfit between the current M1 (x) and reference models was estimated as follows
0,0001
Single charge-discharge cycle [M.Years]
I i
0.04 0.06 0.08
E(dt,x) = \\MC(dt,x) - MREF(x)\\/\\MREF(x)\\,
V die aß, ; VxeX,
(40)
where the least square vector norm is implemented.
These three steps were routinely repeated for different positions of the model vector within the model space. The dynamics of modelling errors vs. single charge-discharge cycle is indicating in Fig. 11.
The stability computer tests show the following.
The combined scheme provides reliable converging of current synthetic data to the reference
Fig. 11. The time slicing error E(dtx) of the combine numerical scheme at the excess pressure (a) and porosity (b) modeling vs. time grid density The modeling error (vertical axes) is given in percent according to formula (40). The vertical bars indicate an aperture of multi-curve averaging (bold curve) built based on 8-model tests.
one when the number of single charging-discharging cycle < 40 Ky per formation time interval burying with maximum sedimentation rate up to 1000m/M.Years. At that, the misfit value between synthetic data computed on such gross grid and the reference prototype is still below 1 % in average. I.e. it 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 the modelling error as a function of the model vector position is not stationary. Depending on the model parameter range, the error estimation can randomly exceed 10 % level15.
The single time cycle threshold for the moderate sedimentation rate ranged between 50-500 m/My is about 100 Ky. Below this level the modelling error (40) is rarely exceed 0.5 % and converging process becomes stationary. Above this level the averaged error gradually decreases (Fig. 11).
Speed
The speed of the combined 3-D solution algorithm is considered here against the classical regular grid 3-D scheme of the relevant forward HD problem solution (Wang, Anderson, 1982). Assume that the same implicit approach based on the ADI method is applied for x, y and z directions at both of the numerical schemes. It allows to estimate the benefit of combining the 1.5-D and 2.5-D solutions.
Let the regular model spatial and time grids be given by
13 For regular grids they are determined by discreteness threshold: means of subdividing on sub-homogeneous intervals, their length and the single grid step inside each of them.
14 1Ky = 1000 years.
15 Such big spikes are associated normally with the discontinuity in burial history: erosion episodes, fault activation interval and so on.
0.2
<
(fxyz = ((xjk): xt = xmin + ihx, i = 0,1,2,..N y = ymin + jhy, j = 0,1,2,..,Ny, zk = khz, k = 0,1,2,...,Nz},
x
f111 + hxNx = xmax, ymin + hyNy = xmax, Kxhz = zmax . Nht = T.
(41)
e?t = (tk: tk = nht, n = 0,1,2,.,.,N - 1,N},
In order to get the 3-D solution at every k-th time slice of the grid cJRt one should get C(«RRx,yz) = NxNy+ NxNz + NyNz solutions of 1-D HD problem along x, y and z directions, respectively.
Let us now get the 3-D combined solution by using the most stable and time consuming mode, where a single time step is equal to a single charging-discharging cycle. Note, that only this particular construction of time grid coincide with the classical approach.
Now one should apply a five-step routine described above at the k-th time slice of the same time grid to get C(®Tx>yueRRz) solutions including:
2MV of 1-D HD problems posed on regular grid oRz (z direction), where 2MV is a total number of 1-D models
required for the target ^area after upscaling of the relevant Earth model; ML(Nx + Ny) of 1-D HD problems posed on regular grids ¿/^ and cJRy for x and y directions, respectively, where ML is a number of 2-D aquifer formations detected within the target Q area after upscaling of the relevant Earth model. Hence, the combined scheme of the 3-D HD problem solution requires for the same target area all together C(®Tx>yueRRz) = MLNx + MLNy + 2MV solutions of the 1-D HD problem on the same time grid. Thus, the reduction coefficient of CPU time rt can be estimated as follows:
Practically, the single dimension of the regular spatial grid cJR has range: n102 - m103 nodes. The upscaling routine allows reducing of grid requirements for combined 3-D modelling in orders (see Chapter 2). Let, for example, the spatial grid requirements of classical 3-D grid approach be the following: Nx = Nz = Nz = 1000. Let for the same target area Q get (after Earth model upscaling) the following requirements for the combined grid: ML = 5 and MV = 50 at the same regular grid aRz for the 1-D problem. Then the reduction coefficient rt according to (42) is equal to 0.00336(6). In other words, the CPU time requirement of combined 3-D numerical scheme constitutes about 0.34 % of the corresponding characteristic of the classical 3-D numerical scheme.
Note, that formula (42) gives an approximate value of reduction, since it does not take into account delay of combined scheme due to interpolating MV of 1-D solutions on ML regular x, y grids {dRx,y} (see STEP 3 of the algorithm above). On the other hand, the potential reserve for additional acceleration of the combined scheme is in extending charging-discharging time cycle.
The combined scheme of 3-D HD modelling is also significantly saving computer memory. Indeed, the memory needs for classical way on regular 3-D grid can be estimated as NxNyNz of real array elements. The same needs for the combined scheme are equal to MLNxNy + MVNz elements. The memory reduction coefficient rM can be estimated as follows:
There are many possible applications of the developed Effective Basin Model approach to the overpressure modelling. Still, there are some of them, where 3-D pressure communication effects are most important for understanding of the observable phenomena and possibility to predict them before drilling. It is clear, that the only 3-D modelling tool is capable to get realistic picture of the present overpressure distribution in subspace. From this point of view it is crucially important to keep 3-D features of the target phenomenon while its computer simulation and be able to run such model in reasonably short single run CPU time. Below we present the results of testing introduced combine 3-D forward modelling operator on practically important centroid effect and pressure compartmentalisation effects16.
6.1. Simulation of centroid effect
The Gulf of Mexico geology settings (Mello, Karner, 1996) was taken as a prototype of Centroid conditions in terms of stratigraphy O lithology links, burial history and present day aquifer geometry.
The basin model was established based on five formation units (Fig. 12) with clay (green), sandstone (brown) and evaporite (salt at the bottom part) lithologies.
16 Readers of the electronic version of this article can also watch dynamics of every modelling example by running video-clip avalible at: http://vestnik.mstu.edu.ru/v08 1 n19/articles/01 sered1.avi; http://vestnik.ms1u.edu.ru/v08_1_n19/articles/01_sered2.avi; http://vestnikmstu.edu.ru/v08_1_n19/articles/01_sered3.avi
rt = C(^x,y,z)/C(©Tx,yU^z) = (MLNx + MLNy + 2Mv)/(Nx Ny + Nx Nz + NyNz). (42)
ru = (MLNxNy + MvNz) / NxNyNz. This gives the value rM = 0.00505 for the same example (Nx = Nz = Nz = 1000 and ML = 5, MV = 50).
(43)
6. Examples
in
-
B Pie Pie
Pia
Pliocene
Zandian RS1
Mes
¡"ortonian
]s sar
e >- ■I
. [M Sev
e emi H Miocene
.anghian
Bur SI
BecmmiK MTTY, mom 8, №1, 2005 2. cmp.5-43
Enin.il historv
2.4
b
4 3. 2 Tirce
1.6 0.8 Burial time [mYears]
c
1.0
1.5
Pressure
2.0 g/c3 [EMW]
d
1G 12.5
Distance [ft] H Ü000)
0
a
Fig. 12. The local basin setting at the test of centroid pressure effect (the Gulf of Mexico example).
Young (a) and old (b) stratigraphy models. Typical 1-D burial model (c) and associated 1.5-D pressure modelling result (d). Cross section (e) through the crest of salt induced dome with synthetic pressure profiles: 1-D model (black line) vs. 3-D young (blue markers) and old (red markers) models. The pseudo well corresponding to c, d fragments is framed.
Fig. 13. The test of centroid effect for young (lower) and old (upper) sedimentation models introduced in Fig. 12a-b. The misfit maps (left) between 1D and combine 3D pressure solutions are given in mPa. The colour excess pressure scales for map window (middle) and section window (right) are common.
A centroid pressure effect occurs at the top part of the steep arc shaped aquifer formation if its average HD conduction capacity and sedimentation time available before the present day allow equalising of excess hydrostatic pressure inside this HD unit (Traugott, 1996). At that, the pressure difference between clay part of the section where lateral pressure transfer is close to zero and good laterally permeable bed could be dramatic (especially at the crest part). Note that the only using of a 3-D pressure communication model can simulate this effect reasonably.
A fast estimation can be made based on simple assumptions.
Indeed, the difference in permeability between the clays and sandy rocks could vary in several orders. According to (Mello, Karner, 1996) the permeability of sandstone rocks at the depth (16-20)-103 ft (5-6.5 km) in the Gulf of Mexico conditions rarely could reach 1milli Darcy (mD). As it was shown (Madatov, Sereda, 2000a,b) the relevant HD conductivity level provides the speed of Darcy flow for pore water movement within the
1.0
—
— 1.2
1.4
1 1.6
1.8
2.0
Pressure in g/cc of EMW
Fig. 14. A typical pore pressure compartmentalisation in the North Sea region. Challenge well trace is highlighted with the pink color
s e
r
a
e r o f e
LO
H
Tertiary M i o ce n e
□ ligo ce ne
E o ce n e
Pale <:• ce n e
Cretaceous Gull K2
Early K1
. ,
M a I m J3
Mesozoic
a Jurassic Dogger J2 :
Lias J1
Tr3
£ ft
2250 2500 2750 3000 3250 3500 3750 4000 4250 4500 4750 5000
d
200 150 100 50 Burial time [mYears]
b
0 1.0 1.5 2.0 g/c3 Pressure [EMW]
c
2.0 3.5
Distance [m] *10000
Fig. 15. The local basin setting at the test of pressure f compartmentalization effects (the North Sea example).
a - stratigraphy models with aquifer marked; b,c - the typical 1-D burial model (b) and associated 1.5-D pressure modelling result (c); d - the target AU subsurface map; e - cross section with offset (black) and
"planed" (pink) wells marked
range of about 10"5-10"3 m per 103 years at the constant pressure gradient 1 Pa/m. Based on 1-D pressure estimations at the flank and crest parts of the dome sample, one can expect 103 Pa/m lateral pressure gradient value as an upper limit during burial history. So, the maximum distance of lateral pore water movement could be estimated for this example as 2.5 km for the young model and 25 km for the old model, respectively. Consequently, the age of aquifer formation on the condition that the other model parameters are fixed can play here a paramount role.
The new 3-D modelling tool can easily check both young and old stratigraphy scenarios for the HD system associated with the dome model in Fig. 12.
Combined pore pressure solution for aquifer formation R7 reveals clear visible centroid effect (Fig. 13, upper) at the old version of the burial history setting. The difference between interpolating of 1-D pressure solutions within this aquifer formation and 3-D combined solutions was checked. The young model case reveals a very little pressure difference between 1-D and 3-D solutions at control points, whereas the old case produces up to 100 % excess pressure difference at flank and top parts.
6.2. Modelling of pressure compartmentalisation
A North Sea basin represents an ideal testing area for investigation and simulating of pressure compartmentalisation phenomena. The mosaic present day pore pressure distribution at the Jurassic reservoir level could be explained by using at least four following reasons (Chiarelli, Duffaud, 1980):
1. Fault-block tectonics on the rift megasequence level can create significant lateral seals between different pressure compartments.
2. Two different types of subsidence during rifting and post-rifting time which leads to quite isolated pressure regimes for upper and deeper Cretaceous unconformity megasequences.
3. Several periods of non-deposition or uplifting and erosion during Jurassic-Early Cretaceous time which could cause - 15600E-05 simultaneous local release of the overpressured pore fluid.
4. Higher temperature gradient in the pre-axial part of the Graben due to different sedimentation and cooling rates. Consequently different start times for the maturation and HC-generation from source rocks and different phase content of the pore fluid observed in different pressure compartments at the present day.
As a result the mechanisms responsible for the pressure development vary from one pressure compartment to another and from the upper level to the deeper level in the section (Fig. 14). More detailed review of the overpressuring mechanisms acting in the region is available in the published literature (Chiarelli, Duffaud, 1980; Buhrig, 1989; Helth et al, 1994).
The commonly good permeable Middle Jurassic oil-gas reservoir reveals the highest sensitivity among Effective Lateral Conduction parameters17 of the formation set included in the area sensitivity analysis. According to the Earth model upscaling strategy (Madatov, Sereda, 2003) it was used as the main aquifer formation responsible for pressure communication to establish combined 3-D hydro dynamics model in the basin scale.
0.19126E-05 i.OOOOOE+OD
c
lip tot
d
0 66285E-03.
Pressure [EftjW units]
e
f
Fig. 16. Results of combine 3-D HD modelling the North Sea example. a - the aquifer HD permeability grid map (in LOG scale) [D]; b - the vector lateral flow grid map; c - the current lateral sink term grid map [Sp.surface]; d - the current vertical sink term grid map [Sp.surface]; e -the current excess hydrostatic pore pressure grid map [g/cc EMW]; f - the current excess hydrostatic pore pressure along section trace (black) in the same scale Click here to get video clip
The lateral sink term.
The grid structure map (see Fig. 15a) of the aquifer top was analysed on presence of discontinuities associated with the fault system. The analysis and further section interpretation reveals several terraces separated by the steeply dipping and vertical faults. The target subsurface of aquifer formation was finally represented on triangular mesh extended with the set of fault nodes reproducing local tectonics features by linear segments of vertical faults (Fig. 15d-e). Some of these faults could be recognised as sealing ones by joint interpretation of subsurface topography and formation pressure data if there are available at the relevant depth interval.
Note, that the quality of lateral seal and consequently the pressure communication at any specific position can only be estimated after careful calibration of relevant fault transmissibility parameters (Madatov, Sereda, 2001). The first step in this direction closely depends on good guess in the seal model parameterisation, the start values of relevant membrane characteristics (Antonellini, Aydin, 1994; Madatov, Sereda, 2001) and appropriate 3-D forward modelling operator.
The 3-D combine modelling results for aquifer unit (marked with arrow in Fig. 15a) are represented in Fig. 16-17.
Fig. 16 represents a work screen of the 3-D modelling program during running on input corresponded to the model setting (Fig. 15). Here, the little white squares on the grids highlight the position of the offset wells where formation models and pore pressure were available. Together with the pseudo wells they form an optimised triangular mesh used for computing of the combined 3-D HD model. The x, y distributions of lateral (Fig. 16c) and vertical (Fig. 16d) sink terms on regular grids aPxy at each time slice of combined excess pressure solution reflect the current balance between the vertical and lateral (Fig. 16b) fluid flows. Note, that it is spatially controlled by lateral conductivity of burying aquifer formation (Fig. 16a). The upward lateral release of the pore water from the deepest compartments had two major directions: to the east flank by multi-steps across fault patterns and opposite, to the west - across the single sealing fault with the big displacement. The east direction provides a main lateral discharge of the overpressured HD system during post activation time interval which is clearly visible in Fig. 15b. As a result, the present day pressure drop across the west fault is proved to be
cn -<
ft 1 '
u Q
LJT:-
tC ■ ■-
G/cc EMW
1 1.5 2.
G/cc EMW
G/cc EMW
G/cc EMW
Distance [m 10000]
a
Fig. 17. Comparison of 1.5-D pressure solutions for the set of wells along section (green line trace) against 3-D combine pressure solution for aquifer formation. a - the section across several pressure compartments. The vertical well paths are black (pink for target well location). The relevant 1.5-D pressure models attached (light red curves). The 3-D pressure solution for aquifer formation is indicating with the blue triangle markers. The difference between 1.5-D and 3-D solution at every single well is highlighted (bold red strip on excess pressure scale); b - the present day grid model of excess pore pressure for aquifer formation with section trace (green) and wells involved into combine 3-D modelling (red circles). The scale of map is given in g/cc of EMW. The detected lateral seals are highlighted.
1 1.5
2
1 1.5
2
1 1.5
2
sharpest. Some non-balanced lateral drainage in the east direction resulted with residual excess pressure detected on the relevant flanks. Due to steep dipping and old age of aquifer formation the appearance of centroid effects at the highest tectonically sealed traps are possible.
Finally, the disagreement between the 1.5-D and combined 3-D pore pressure modelling results was analysed (Fig. 17). The maximum excess pressure mismatch was detected for the deepest compartment. At the assumed target well location it reached about 100 % (Fig. 17a). Clearly, such severe mistake in pore pressure pre-drill estimations is due to restricted validity of 1.5-D approach for the areas where 3-D pressure communication effects play a paramount role during the burial history.
7. Conclusion
The review and analysis of the existing BM approaches to the problem of overpressure simulation reveals a big potential for creation of optimally simple tool 3-D overpressure modelling oriented on the further calibration against real data and use for prediction. Briefly the achieved results can be reduced to the following two statements:
• The existing level of computer geo-science technologies allows the possibility to produce 3-D basin models based on rather dense grid and detailed specification. In connection with pore pressure prediction purposes such kind of computer models are unnecessarily heavy and hard to use. This is especially true for the real time updating of pre-drill pore pressure prediction during drilling.
• The developed combined approach to 3-D HD modelling in basin time scale reduces requirement in computer speed and memory in orders without loosing data constrained modelling accuracy. It implies introducing unique for the area lithology and tectonics specifications versus stratigraphy column for the area. It allows compact attribution of the Earth model via four independent formation parameters and a set of the fault transmissibility characteristics. All of them, including the total number of formation units in area formation column and fault segments accounted in the present day tectonics picture, are tuneable during calibration. At that, the aquitard lithology units for tectonically relaxed basin conditions18 are considered as essentially 1-D parts of the sedimentary rock section in terms of basin scale hydrodynamics. In contrast, the aquifer formation units effectively work as lateral pressure connectors during the whole sedimentation history. A combined 3-D HD modelling operator developed as a part of general Earth model upscaling strategy (Madatov, Sereda, 2003) is capable to account for the main 3-D pressure communication phenomena (centroid, dipping aquifer beds, pressure drop across the sealing faults, etc.).
References
Alberty M.W., McLean R.M. Emerging trends in pressure prediction. Paper presented on Offshore Technology
Conference held in Houston, Texas, USA, May 5-8, 2003. Allen P.A., Allen J.R. Basin analysis principles and applications. Blackwell scientific publication, Oxford, London, p.393, 1990.
Antonellini M., Aydin A. Effect of faulting on fluid flow in porous sandstones: Petrophysical properties. The
American Assoc. of Petroleum Geologists Bull., v.78, N 3, p.355-377, 1994. Antonellini M., Aydin A. Effect of faulting on fluid flow in porous sandstones: Geometry and spatial
distribution. The American Assoc. of Petroleum Geologists Bull., v.79, N 5, p.642-671, 1995. Athy L.F. Density, porosity and compaction of sedimentary rocks. The American Assoc. of Petroleum Geologist Bull., v.14, p.1-24, 1930.
Audelt D.M., McConnell J.D.C. Forward modelling of porosity and pore pressure evolution in sedimentary
basins. Basin Research, N 4, p.147-162, 1994. Baldwin B., Butler C.O. Compaction curves. The American Assoc. of Petroleum Geologists Bull., v.69, N 4, p.622-626, 1985.
Barker C. Calculated volume and pressure changes during the thermal cracking of oil into gas in reservoirs. The
American Assoc. of Petroleum Geologists Bull., v.74, p.1254-1261, 1990. Barker C. Generation of anomalous internal pressures in source rocks and its role in driving petroleum
migration. Revue De L'Institut Francais Du Petrole, v.43, p.3, 1988. Barker C. Thermal modelling of petroleum generation: Application and theory. Elsevier, Amsterdam - New
York - Oxford - Shannon - Tokyo, 395 p., 1996. Bear J., Bach mat Y. Introduction to modelling of transport phenomena in porous media. Kluwer Academic Publishers, London, p.553, 1991.
18 There is no any significant horizontal stress activity during burial history.
Berdnikov N.I., Golikova G.V., Grossheim A.V. Application of effective seismic model concept. M., Nedra, 184 p., 1985 (in Russian).
Biot M.A. Theory of propagation of elastic waves in a fluidsaturated porous solid. I. Low-frequency range. Appl. Phys, v.13, p.35-40, 1956.
Bowers G.L. Pore pressure estimation from velocity data: Accounting for overpressure mechanisms besides
undercompaction. SPE, Drilling & Completion, N 6, p.88-95, 1995. Bredehoeft J.D., Djevanshir R.D., Belitz K.R. Lateral fluid flow in compacting sand-shale sequence: South
Caspian Basin. The American Assoc. of Petroleum Geologists Bull., v.72, N 4, p.416-424, 1998. Buhrig C. Geopressured Jurassic reservoirs in the Viking Graben: Modelling and geological significance.
Marine and Petroleum Geology, v.6, p.35-48, 1989. Buland A., Omre H. Bayesian seismic inversion and estimation in a spatial setting. Paper D-39 presented on
EAEG conference, Stavanger, June 4-7, 2003. Campbell A., Elam S., Lahann R., Patmore S. Pressure prediction from seismic derived velocities. Leacon learned at 204/19-1, West of Shetlands, UK. Paper presented on "Overpressure 2000" workshop, 4-6 April, London, 2000.
Chapman D.S., Keho T.H., Bauer M.S., Picard M.D. Heat flow in the Uinta basin determined from bottom hole temperature data. Geophysics, v.49, N 4, p.453-466, 1984.
Charlez P. Rock mechanics. Theoretical Fundamentals, Editions Thecnip., Paris, v.1, 333 p., 1991. 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. Claerbout J.F. Imaging the Earth's interior. Blackwell Scientific Publishing. Oxford, 407 p., 1985. Dobrynin V.M., Serebryakov V.A. Methods of prediction of abnormally high pore pressure. Moscow, Nedra, 232 p., 1978 (in Russian).
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. Domenico P.A., Palciauskas V.V. Thermal expansion of fluids and fracture initiation in compacting sediments.
Geol. Soc. Am. Bull, v.90, p.953-979, 1979. Dutta N.C. Shale compaction, burial diagenesis and geopressure: Dynamic model, solution and some results. In "Thermal modelling in sedimentary basins" collection of papers from 1-st Exploration Research Conference. IFP, France, p.173-184, 1985. Eaton B.A. The equation for geopressure prediction from well logs. Paper No. 5544 presented on SPE Conference, 1975.
England W.A., Mackenzie A.S., Mann D.M., Quigley T.M. The movement and entrapment of petroleum
fluids in subsurface. Journal of the Geological Society, London, v.144, p.120-139, 1987. Fault mechanics and transport properties of rock. Edited by B. Evans, Academic Press, London, 523 p., 1992. Fox D.R. Estimation and inference in inverse methods for kinetic models of hidrocarbon generation.
Mathematical Geology, v.27, N 8, p.923-938, 1995. Gangi A.F., Berg R.R. Primary migration by oil generation fracturing in law permeable source rocks. The
American Assoc. of Petroleum Geologists Bull., v.81, N 3, p.424-443, 1997. Goffe W.L., Ferrier G.D., Rodgers J. Global optimization of statistical functions with simulated annealing.
J. Econometrics, v.60, p.65-100, 1994. Goldberg. Genetic algorithms in search optimization and machine learning. Addison Wesley Publ. Co, p.281-292, 1989.
Guidish T.M., Kendall C.G., Lerch I., Toth D.J., Yarzab R.F. Basin evaluation using burial history calculations:
An overview. The American Assoc. of Petroleum Geologists Bull., v.69, N 1, p.92-105, 1985. Gurevitch A.E., Kraichik M.S., Batygina N.B. A geo-fluids' pressure. Leningrad, NEDRA, p.223, 1987 (in Russian).
Helth A.E., Walsh J.J., Watterson J. Estimation of effects of sub-seismic sealing faults on effective
permeabilities in sandstone reservoirs. Norwegian Institute of Technology, p.327-347, 1994. Huang Z., Williamson M.A., Katsube J. Permeability prediction with artificial neutral networks modelling in
the Venture gas field, offshore eastern Canada. Geophysics, v.61, N 2, p.422-435, 1996. Hunt J.M. Generation and migration of petroleum from abnormally pressured field compartments. The
American Assoc. of Petroleum Geologists Bull., v.74, p.1-12, 1990. Ingberg L. Very fast simulated re-annealing. Math., Comp. Mod., v.12, p.967-973, 1989. Jouset T.P. Migration of hydrocarbons through the sedimentary rocks. M., Nedra, p.188, 1986 (in Russian).
Katahara K. Analysis of overpressure on Gulf of Mexico Shelf. Paper presented on Offshore Technology Conference held in Houston, Texas, USA, May 5-8, OTC 15293, 2003.
Katsube N., Caroll M.M. The role of Terzaghi's effective stress in linearly elastic deformation. Journal of Energy Resources Technology, v.105, p.509-511, 1983.
Katz S.A., Chilingar G.V., Gurevich A., Khilyuk L. Bi-linear models for simultaneous estimation of formation pressure and lithological characteristics in interbedded sands and shales. J. of Petroleum Science and Engineering, v.12, p.37-48, 1994.
Knott S.D. Fault seal analysis in the North Sea. The American Assoc. of Petroleum Geologists Bull., v.77, N 5, p.778-792, 1993.
Landro M. Discrimination between pressure and fluid saturation changes from time-laps seismic data. Geophysics, v.66, p.836-844, 2001.
Leonard R.C. Distribution of sub-surface pressure in the Norwegian Central Graben and application for exploration. From Petroleum Geology of North-West Europe of the 4-th Conference, p.1295-1303, 1993.
Lerch I. Inversion of dynamical indicators in quantitative basin analysis models. 1. Theoretical considerations. Mathematical Geology, v.23, N 6, p.817-832, 1991.
Lerch I. Theoretical aspects of problems in basin modelling. In "Basin Modelling Advances and Applications", Norwegian Petroleum Society, Special publication, Elsevier, Amsterdam, p.35-65, 1990.
Liu G., Roaldset E. A new decompaction model and its application to the northern North Sea. First break, v. 12, N 2, February, p.81, 1994.
Luo X., Vasseur G. Contributions of compaction and aquathermal pressuring to geopressure and the influence of environmental conditions. The American Assoc. of Petroleum Geologists Bull., v.76, N 10, p.1550-1559, 1992.
Madatov A.G., Sereda A.-V.I. Computing aspects of the solution of some inverse problems in geophysical investigations. Abstract of papers to be presented at the International Conference "Northern Universities". MSTU, Murmansk, p.57-58, 1997.
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 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., Doyle E.F., Helle H.B. The "1.5-D" inversion approach to the pore pressure evaluation. Concept and application. Paper presented on Workshop "Compaction and Overpressure Current Research", 9-10 December 1996, Institute Français du Petrole, Paris, France, 1996.
Madatov A.G., Sereda V.-A.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 V.-A.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 American Assoc. of Drilling Engineers (AADE) Forum "Pressure regimes in sedimentary basins and their prediction", 2-4 September 1998, Houston, Texas, USA, 1999.
Madatov A.G., Sereda V.-A.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.
Mallick S. Some practical aspects of prestack waveform inversion using a genetic algorithm: An example from the east Texas Woodbine gas sand. Geophysics, v.64, N 2, p.326-336, 1999.
Mann D.M., Mackenzie A.S. Prediction of pore fluid pressure in sedimentary basins. Marine and Petroleum Geology, v.7, N 2, p.55-65, 1990.
Manzocchi T., Walsh J.J., Nell P., Yielding G. Fault transmissibility multipliers for flow simulation models. Fault Analysis Group. The University of Liverpool Department of Earth Sciences, August, p. 1-37, 1998.
McCain W. The properties of petroleum fluids. Tulsa, Oklahoma, 548 p., 1990.
Mello U.T., Karner G.D. Development of sediment overpressure and its effect on thermal maturation: Application to the Gulf of Mexico basin. The American Assoc. of Petroleum Geologists Bull., v.80, N 9, p.1367-1395, 1996.
Miall A.D. Principles of sedimentary basin analysis. Springer-Verlag, Berlin, p.616, 2000.
Migration of hydrocarbons in sedimentary basins. The collection of paper from 2-nd Exploration Research Conference. Edited by Doligez B., IFP, France, p.681, 1987.
Miller T.W. New insights on natural hydraulic fractures induced by abnormally high pore pressure. The American Assoc. of Petroleum Geologists Bull., v.79, N 7, p.1005-1018, 1995.
Mitra S., Beard W.C. Theoretical models of porosity reduction by pressure solution for well-sorted sandstones. J. Sediment. Petrol, v.50, p.1347-1360, 1980.
Moore C.H. Carbonate diagenesis and porosity. Elsevier, Amsterdam, 338 p., 1989.
Morley Ch.K. Patterns of displacement along large normal faults: Implications for basin evolution and fault propagation, based on example from East Africa. The American Assoc. of Petroleum Geologists Bull., v.83, N 4, p.613-634, 1999.
Mouchet J.P., Mitchell A. Abnormal pressures while drilling, Manuels techniques. Elf Aquitaine, Boussens, 286 p., 1989.
Mudford B.S. A one-dimensional, two-phase model of overpressure generation in the Venture gas field, offshore Nova Scotia. Bull. of Canadian Petroleum Geology, v.38, N 2, p.246-258, 1990.
Oil and Gas Geology handbook. M., Nedra, 480 p., 1984 (in Russian).
Ortoleva P.J. Basin compartments and seals. The American Assoc. of Petroleum Geologists, Memoir 61, Tulsa, Oklahoma, 1999.
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.
Patterson. Artificial neural networks: Theory and applications. Prentice Hall, 477 p., 1996.
Reike H.H., Chilingarian G.V. Compaction of argillaceous sediments. London, New York, 424 p., 1974.
Rittenhouse G. Pore space reduction by solution and segmentation. The American Assoc. of Petroleum Geologists Bull., v.55, p.80-91, 1971.
Rowan M.G. A systematic technique for the sequential of salt structures. Tectonophysics, v.228, p.369-381, 1993.
Rybach L., Buntebarth G. The variation of heat generation, density and seismic velocity with rock type in the continental lithosphere. Tectonophysics, v.103, p.335-344, 1984.
Samarsky A.A., Gulin A.V., Numerical methods. M., Nauka, 430 p., 1989 (in Russian).
Schneider F., Bouteca M., Vasseur G. Validity of the porosity / effective-stress concept in sedimentary basin modelling. First break, v.12, N 6, p.321-326, 1994.
Schneider F., Potdevin J.L., Wolf S., Faille L. Modele de compaction elastoplastique et viscoplastique pour simulateur de bassins sedimentaires. IPFRevue, v.49, N 2, p.141-148, 1994.
Seismic Stratigraphy. Edited by Ch.E. Payton. Memoir of the American Assoc. of Petroleum Geologists, Tulsa, Oklahoma, p.839, 1977.
Sen M.K., Stoffa P.I. Bayesian inference. Gibbs' sampler and uncertainty estimation in geophysical inversion. Geoph. Prospecting, v.44, p.313-350, 1996.
Shi Y.L., Wang C.Y. Pore pressure generation in sedimentary basins: Overloading versus aquathermal. J. Geophys. Res, v.91, p.2150-2162, 1986.
Smith J.E. The dynamics of shale compaction and evolution of pore-fluid pressures. Math. Geol., v.3, p.239-263, 1973.
Spencer C.W. Hydrocarbon generation as a mechanism for overpressuring in Rocky-Mountain region. The American Assoc. of Petroleum Geologists Bull., v.71, p.368-388, 1987.
Spiers C.J., Schutjens P.M., Brzesowsky R.H., Peach C.J., Liezenberg J.L., Zwart H.J. Experimental determination of constitutive parameters governing creep of rocksalt by pressure solution. In: R.J. Knipe and E.H. Rutter (Editors), Deformation Mechanisms, Rheology and Tectonics. Geol. Soc. London, Spec. Publ, v.54, p.215-227, 1990.
Straughan B., Greve R., Ehrentraut H., Wang Y. Continuum mechanics and applications in geophysics and the environment. Springer-Verlag, Berlin Heidelberg, p.393, 2001.
Stump B.B., Flemings P.B., Finkbeiner T., Zoback M.D. Pressure difference between overpressured sands and bounding shales of the Eugene island 330 Field (Offshore Louisiana, USA) with implication for fluid flow induced by sediment loading. Paper presented on Workshop "Overpressure in Petroleum Exploration", Pau, April, 1998.
Stuve K. Geodynamics of the lithosphere (An introduction). Springer-Verlag, Berlin Heidelberg, p.449, 2002.
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.
Swarbrick R.E., Osborne M.J. The nature and diversity of pressure transition zones. Petroleum Geoscience, v.2, p.111-116, 1996.
Sweeney J.J., Burnham A.K. Evaluation of a simple model of vitrinite reflectance based on chemical kinetics.
The American Assoc. of Petroleum Geologist Bull., v.74, N 10, p.1559-1570, 1990.
Sweeney J.J., Burnham A.K., Braun R.L. A model of hydrocarbon generation from type I Kerogen: Application to Uinta Basin, Utah. The American Assoc. of Petroleum Geologist Bull., v.71, p.967-985, 1987.
Tarantola A. Inverse problem theory: Methods for data fitting and model parameter estimation. Elsevier, Netherlands, p.386, 1987.
Terzaghi K., Peck R.B. Soil mechanics in engineering practice. Wiley, New York, 566 p., 1948.
Thermal modelling in sedimentary basins. The collection of paper from 1-st Exploration Research Conference. Edited by J.Burruse, IFP, France, p.600, 1985.
Thorne J.A., Watts A.B. Quantitative analysis of North Sea subsidence. The American Assoc. of Petroleum Geologists Bull., v.73, N 1, p.88-116, 1989.
Tissot B.P., Welte D.H. Petroleum formation and occurrence. Springer-Verlag, New York, p.538, 1978.
Traugott M. The pore pressure centroid concept: Reducing drilling risk. Compaction and Overpressure Current Research, Abst., 9-10, December, IFP, Paris, 1996.
Ungerer P. Modelling of petroleum generation and expulsion - an update to recent reviews. Basin Modelling: Advances and Applications, NPF Special Publ., v.3, p.219-232, 1993.
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.
Verweij J.M. Hydrocarbon migration systems analysis. Development in Petroleum Science, v.35, p.276, 1993.
Wang H.F., Anderson M.P. Introduction to groundwater modelling. Finite difference and finite element methods. Academic Press Inc., 237 p., 1982.
Waples D.W., Kamata H. Modelling porosity reduction as a series of chemical and physical processes. NPF Elsevier, Amsterdam, Special Publication, v.3, p.303-320, 1993.
Weedman S.D., Brantley S.L., Shiraki R., Poulson S.R. Diagenesis, compaction and fluid chemistry modeling of sandstone near a pressure seal: Lower Tuscaloosa formation, Gulf Coast. The American Assoc. of Petroleum Geologists Bull., v.80, N 7, p.1045-1084, 1998.
Yardley G. Can lateral transfer explain the high pressure in the Central North Sea? Paper presented on Workshop "Overpressure in Petroleum Exploration", Pau, April, 1998.
Yassir N.A., Bell J.S. Relationships between pore pressure, stress and present day geodynamics in the Scotian shelf, Offshore Eastern Canada. The American Assoc. of Petroleum Geologists Bull., v.78, N 12, p.1863-1880, 1994.
Yu Z., Lerch I. Modelling abnormal pressure development in sandstone/shale basins. Marine and Petroleum Geology, v.13, N 2, p.179-193, 1996.
Yu Z., Lerch I., Bour Q. Inversion of dynamical indicators in quantitative basin analysis models. II Synthetic multiwell information and two-dimensional case history. Mathematical Geology, v.27, N 1, p.41-68, 1995.
Zhao K., Lerch I. Inversion of dynamical indicators in quantitative basin analysis models. II Synthetic test and case history. Mathematical Geology, v.25, N 2, p.107-120, 1993.
