Computer implemented method for generating a subsurface rock and/or fluid model of a determined domain

Abstract

A computer implemented method is provided for generating a subsurface rock and/or fluid model on a determined domain. The subsurface rock and/or fluid model generated by the invention combines acoustic data and data obtained by an electromagnetic survey in such a way that the resulting model avoids the need of generating an image of the resistivity with the low resolution imposed by the known techniques inverting EM electromagnetic survey data.

Claims

1. A computer implemented method for analysing a seismic survey dataset and an electromagnetic survey dataset to generate a sub-surface rock and/or fluid model of a specified domain, said method comprising: a) generating an acoustic and/or elastic model of the sub-surface domain based on a seismic survey data; b) generating a rock and/or fluid model populated with rock and/or fluid properties based on the elastic model compatible with acoustic data or on the acoustic model and, wherein the rock and/or fluid model is at a resolution of the acoustic model or the elastic model; c) iteratively, carrying out the following steps: determining a misfit according to: i. determining the acoustic or elastic response of the rock and/or fluid model and, determining a first misfit between said acoustic or elastic response and the acoustic and/or elastic model; or ii. determining the electric properties of the rock and/or fluid model wherein the electric properties are at the scale of the acoustic or elastic model and, determining the electromagnetic response of the rock and/or fluid model according to the electric properties under a simulation of electromagnetic response corresponding to the electromagnetic survey dataset, the simulation being a numerical simulation of the electric properties derived from the rock and/or fluid model for determining the electromagnetic response according to the rock and/or fluid properties of said model and to an electromagnetic source reproducing the electromagnetic fields of the electromagnetic survey and, determining a second misfit between said electromagnetic response and the electromagnetic survey dataset at the locations of a set of electromagnetic sensors of the electromagnetic survey; or iii. determining both the first misfit and the second misfit according to step i) and step ii), respectively; updating the properties of the rock and/or fluid model; wherein the iterations are carried out until the first misfit, if determined, is under a first predetermined value and the second misfit, if determined, is under a second predetermined value.

2. The method according to claim 1, wherein determining the acoustic or elastic response of the rock and/or fluid model of step i) comprises a numerical simulation of the rock and fluid model for determining the acoustic and/or elastic response according to the rock and/or fluid properties of said model.

3. The method according to claim 1, wherein the rock properties include at least one of porosity, clay content, and water saturation.

4. The method according to claim 1, wherein elastic properties of the elastic model include at least one of a P-impedance, a S-impedance, density, Poisson's ratio, Lamé impedances lambda-rho (λρ) and mu-rho (μρ), and Vp/Vs ratio, wherein Vp is P-wave velocity and Vs is S-wave velocity.

5. The method according to claim 1, wherein the elastic response is: generated by using a stiff sand model; or generated by using a soft sand model; or generated by using a rock physics model linking elastic properties, rock properties, and fluid properties; or generated statistically from a relationship between elastic/acoustic response and rock and fluid properties measured in well log data; or generated by any combination thereof.

6. The method according to claim 1, wherein the electric properties include at least one of horizontal resistivity, vertical resistivity, bed parallel resistivity, bed perpendicular resistivity, and anisotropy ratio.

7. The method according to claim 1, wherein the elastic response is: generated by using Archie's model; or generated by using the Simandoux model; or generated by using a rock physics model linking electric properties, rock properties, and fluid properties; or generated statistically from a relationship between electric response and rock and fluid properties observed in well log data; or generated by any combination thereof.

8. The method according to claim 1, wherein an iterative step d) is according to sub-step iii) and wherein the properties of the rock and/or fluid model being updated are at least one of the rock and/or fluid properties.

9. The method according to claim 1, wherein the rock and/or fluid model is a rock and fluid model comprising rock properties and fluid properties and, wherein an iterative step d) combines a sequential process with a first group of iterations according to sub-step i) updating at least one of the properties of the rock and fluid model to minimize the first misfit and a second group of iterations according to sub-step ii) updating at least one of the fluid properties of the rock and fluid model to minimize the second misfit.

10. The method according to claim 1, wherein the rock and/or fluid model is a rock and fluid model comprising rock properties and fluid properties and, wherein the rock and fluid model comprises a plurality N.sub.e of individual rock and fluid models m.sub.j, j=1 . . . N.sub.e each with different rock and fluid property attributes, being each individual rock and fluid model m.sub.j at the resolution of the acoustic model or the elastic model and, being each rock and fluid property attribute dependent on the elastic model or the acoustic model; updating the properties of each of the individual rock and fluid models m.sub.j using a Bayesian updating process such that the first misfit, the second misfit or both are minimised; and, step d) comprises: calculating the average of the plurality of individual rock and fluid models m.sub.j, j=1 . . . N.sub.e, calculating uncertainty over the plurality of individual models, preferably the standard deviation, to provide an estimate of model uncertainty; and making available both values in the model ensemble, the average model and the model uncertainty.

11. The method according to claim 10, wherein the updating process perturbs the properties of the rock and fluid model comprising a predetermined N.sub.a number of updating sub steps, where the amplitude of the updates is decreased from the first updating substep to the last updating substep.

12. The method according to claim 11, wherein the updates in the individual rock and fluid model m.sub.j, is performed using the following relationship:
m.sub.j.sup.p:=m.sub.j.sup.f+C.sub.MD.sup.f(C.sub.DD.sup.f+α.sub.iC.sup.d).sup.−1(d.sub.uc,j−d.sub.j.sup.f) wherein “:=” denotes the updating process as an update from former values, the superscript “p” denotes the updated value of the property, the superscript “f” denotes the previous model value of the property, C.sub.MD.sup.f is the cross covariance matrix between the prior vector of model parameters m.sup.f and the vector of predicted data d.sup.f, C.sub.DD.sup.f is the autocovariance matrix of predicted data, α.sub.i a number of updates being α.sub.i∈[1, N.sub.a] and, d.sub.f.sup.j is the predicted response from the m.sub.j individual rock and fluid model, and d.sub.uc is given by
d.sub.uc=d.sub.obs+√{square root over (α.sub.i)}C.sub.D.sup.1/2z.sub.d where z.sub.d˜N(0,I.sub.N.sub.d) wherein d.sub.uc is the updated vector of observed data, d.sub.obs the observed property, that is, the elastic/acoustic properties determined from the seismic survey data and the electromagnetic data observed in the CSEM electromagnetic survey data, α.sub.i a number of perturbing updates being α.sub.i∈[1, N.sub.a] and, C.sub.D the covariance matrix, z.sub.d a sampled value of a Gaussian distribution.

13. The method according to claim 10, wherein the plurality of individual rock and fluid models is generated as follows: step d) comprises generating a plurality of individual rock and fluid models by stochastic sampling of the reservoir properties dependent only on the acoustic or elastic model; steps c.i-c.iii are carried out for the plurality of individual rock and fluid models and wherein step d) further comprises selecting the subset of the generated individual rock and fluid models for which the electromagnetic response under a simulation of the EM electromagnetic survey match the electromagnetic survey data within a user defined tolerance.

14. The method according to claim 10, wherein the iteration over the plurality of individual rock and fluid models is computed in a parallel manner by grouping the computation of one or more individual rock and fluid models in a separated processor.

15. A computer program product comprising instructions which, when the program is executed by a computer, causes the computer to carry out a method according to claim 1.

Description

DESCRIPTION OF THE DRAWINGS

(1) These and other features and advantages of the invention will be seen more clearly from the following detailed description of a preferred embodiment provided only by way of illustrative and non-limiting example in reference to the attached drawings.

(2) FIG. 1 This figure shows a data processing system for carrying out a method according to the invention.

(3) FIG. 2 This figure shows a schematic representation of the method according to a first embodiment of the invention.

(4) FIG. 3 This figure shows a second embodiment of the invention wherein the elastic iteration and the electric iteration is carried out separately.

(5) FIG. 4a These figures show a plurality of rock and fluid properties of a simple 1D rock and fluid model as an example of the second embodiment.

(6) FIG. 4b These figures show a plurality of acoustic properties (Acoustic impedance, Ip) and elastic properties (Poisson's ratio, PR) in a simple 1D model as an example of the second embodiment.

(7) FIGS. 5a-5d These figures show the evolution of the porosity, clay content and water saturation of the 1D example of FIGS. 4a and 4b according to the iterative process according to the second embodiment.

(8) FIG. 6 This figure shows the evolution of the first misfit (or elastic misfit) and the second misfit (or electric misfit) for the sequence of iterations.

(9) FIGS. 7a,7b These figures show a schematic of a third embodiment of the invention wherein the joint inversion of the rock and fluid properties (lithology of the rock and facies) and the electromagnetic survey data is undertaken using an ensemble method representing statistically the rock and fluid model by means of a plurality of individual rock and fluid models.

(10) FIG. 8a This figure shows a synthetic offshore CSEM amplitude and phase dataset (as an example) generated by forward modelling from a well log based resistivity model. Hereinafter, any synthetic data will be treated as being survey data collected from a seismic survey, a CSEM survey or from both.

(11) FIG. 8b This figure shows the P-impedance and the Poisson ratio generated as synthetic elastic properties from the same well as the CSEM electromagnetic data.

(12) FIG. 9a, 9b These figures show the results of the method according to the third embodiment.

(13) FIG. 10 This figure shows a schematic representation of the steps of a forth embodiment of the method of the invention.

(14) FIG. 11a, 11b These figures show the porosity, water saturation and clay content determined by a rock and fluid model comprising a large amount of individual rock and fluid models sampled from seismic data and a selection of said individual rock and fluid models showing a response being compatible with the synthetic CSEM electromagnetic data resulting in the lowest second misfit.

(15) FIG. 12a This figure shows the phase and amplitude response of the rock and fluid model according to the fourth embodiment.

(16) FIG. 12b This figure shows the misfit between the observed CSEM response and the response calculates from 100 model realisations.

DETAILED DESCRIPTION OF THE INVENTION

(17) As will be appreciated by one skilled in the art, aspects of the present invention may be embodied as a system, a method or a computer program product. Accordingly, aspects of the present invention may take the form of an entirely hardware embodiment, an entirely software embodiment (including firmware, resident software, micro-code, etc.) or an embodiment combining software and hardware aspects that may all generally be referred to herein as a “circuit,” “module” or “system.” Furthermore, aspects of the present invention may take the form of a computer program product embodied in one or more computer readable medium(s) having computer readable program code embodied thereon.

(18) Any combination of one or more computer readable medium(s) may be utilized. The computer readable medium may be a computer readable signal medium or a computer readable storage medium. A computer readable storage medium may be, for example, but not limited to, an electronic, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any suitable combination of the foregoing. More specific examples (a non-exhaustive list) of the computer readable storage medium would include the following: an electrical connection having one or more wires, a portable computer diskette, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or Flash memory), an optical fiber, a portable compact disc read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination of the foregoing. In the context of this document, a computer readable storage medium may be any tangible medium that can contain, or store a program for use by or in connection with an instruction execution system, apparatus, or device.

(19) A computer readable signal medium may include a propagated data signal with computer readable program code embodied therein, for example, in baseband or as part of a carrier wave. Such a propagated signal may take any of a variety of forms, including, but not limited to, electro-magnetic, optical, or any suitable combination thereof. A computer readable signal medium may be any computer readable medium that is not a computer readable storage medium and that can communicate, propagate, or transport a program for use by or in connection with an instruction execution system, apparatus, or device.

(20) Program code embodied on a computer readable medium may be transmitted using any appropriate medium, including but not limited to wireless, wireline, optical fiber cable, RF, etc., or any suitable combination of the foregoing.

(21) Computer program code for carrying out operations for aspects of the present invention may be written in any combination of one or more programming languages, including an object oriented programming language such as Java, Smalltalk, C++ or the like and conventional procedural programming languages, such as the “C” programming language, Fortran, or similar programming languages. The program code may execute entirely on the user's computer, partly on the user's computer, as a stand-alone software package, partly on the user's computer and partly on a remote computer or entirely on the remote computer or server. In the latter scenario, the remote computer may be connected to the user's computer through any type of network, including a local area network (LAN) or a wide area network (WAN), or the connection may be made to an external computer (for example, through the Internet using an Internet Service Provider).

(22) Aspects of the present invention are described below with reference to illustrations and/or diagrams of methods, apparatus (systems) and computer program products according to embodiments of the invention. It will be understood that each illustration can be implemented by computer program instructions. These computer program instructions may be provided to a processor of a general purpose computer, special purpose computer, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create means for implementing the functions/acts specified in the flowchart and/or block diagram block or blocks.

(23) These computer program instructions may also be stored in a computer readable medium that can direct a computer, other programmable data processing apparatus, or other devices to function in a particular manner, such that the instructions stored in the computer readable medium produce an article of manufacture including instructions which implement the function/act specified in the flowchart and/or block diagram block or blocks.

(24) The computer program instructions may also be loaded onto a computer, other programmable data processing apparatus, or other devices to cause a series of operational steps to be performed on the computer, other programmable apparatus or other devices to produce a computer implemented process such that the instructions which execute on the computer or other programmable apparatus provide processes for implementing the functions/acts specified in the flowchart and/or block diagram block or blocks.

(25) Turning now to the drawings and more particularly, FIG. 1 shows an example of a system 100 for generating a subsurface rock and/or fluid model of a determined domain starting from acoustic field data (AD) comprising a model of acoustic attributes derived from seismic survey data, and EM electromagnetic data, both data collected from separate surveys in the same geographical area, according to a preferred embodiment of the present invention. The condition of being the two data collected from separate surveys in the same geographical area applies to all the embodiments.

(26) A preferred computing system 100 includes one or more computers 102, 104, 106 (3 in this example), coupled together, e.g., wired or wirelessly over a network 108. The network 108 may be, for example, a local area network (LAN), the Internet, an intranet or a combination thereof. Typically, the computers 102, 104, 106 include one or more processors, e.g., central processing unit (CPU) 110, memory 112, local storage 114 and some form of input/output device 116 providing a user interface. The local storage 114 may store the acoustic data and the EM electromagnetic data being accessible by the plurality of computers 102, 104, 106, processing in parallel a plurality of individual rock and fluid models, each individual process being processed by each computer 102, 104, 106, in such a manner that the combination of the plurality of individual rock and fluid models provides an estimation of the uncertainty as it will be disclosed in an specific embodiment.

(27) The present invention provides a method for generating a subsurface rock and/or fluid model of a determined domain. In this embodiment the domain is located under the surface of the seabed and the method starts from data retrieved from two surveys: a first seismic survey carried out by using vibrators, airguns or another acoustic source and acoustic sensors such as geophones or hydrophones distributed over the seabed or in the water column above it and, a second CSEM electromagnetic survey carried out by using a plurality of electromagnetic sensors, in this specific embodiment distributed according to a regular grid, over the surface of a seabed for reading the electromagnetic field generated by a horizontal electric dipole antenna that is flown over the seabed at a height of in general 30-50 m; wherein the electromagnetic field deployed by the antenna is chosen based on the estimated rock and fluid structure located below the seabed, and is recorded by the seafloor sensors.

(28) FIG. 2 gives an overview of the proposed method. Serial numbers in brackets will identify the steps represented in figures with the same number in a circle; that is, these number are not reference numbers to elements, they just indicate the order of execution of the steps of the method.

(29) Seismic data collected from the seismic survey are first inverted using one of the many well-known methods available (1). Preferably pre-stack data are used so that elastic attributes (P-impedance and S-impedance and derived attributes) are calculated. Post-stack data could be used for the inversion, yielding only acoustic properties or attributes. In this embodiment, the data are then converted from time to depth, again using well known algorithms (2). The result is an elastic model comprising a model with elastic properties or attributes, defined in depth, at seismic resolution, to be used as input to the proposed method (3).

(30) The second input to the proposed workflow in this embodiment is processed frequency domain CSEM electromagnetic data collected from the CSEM electromagnetic survey, comprising, in this embodiment, amplitude and phase of the electric and (in some cases) magnetic field, as a function of source-receiver separation, geometry and signal frequency (4). This could be obtained from deep-towed nodal CSEM methods, towed streamer CSEM methods or one of the other CSEM acquisition methods that are well known. Although in this embodiment, frequency domain data are employed, the method could equally be applied to time domain data, in which case the CSEM survey data would comprise decay curves of the electric or magnetic field with time measured at each of the sensors. Similarly the method could be applied to natural source data collected offshore, or to natural source or CSEM data collected onshore or in a borehole.

(31) The next input to the proposed iterative method is a starting rock and/or fluid model, in this embodiment comprising lithology and fluid properties (for example porosity, clay content and water saturation, S.sub.w) and seismically defined facies (5).

(32) The rock and/or fluid model could be a 1D depth profile, a 2D section or a 3D volume. This rock and/or fluid model can be derived from seismic quantitative interpretation processes such as the multi-attribute rotation scheme (MARS), statistical rock physics approaches, or petrophysical inversion of elastic properties only, or other known approaches for seismic quantitative interpretation.

(33) The core of the iterative method is depicted in the dashed box and labelled as (6) in FIG. 1, which represents the joint inversion of seismically derived elastic attributes, and CSEM amplitude and phase data according to the invention. A number of embodiments of this inversion will also be disclosed.

(34) In a first embodiment, shown in FIG. 2, the rock and fluid model (7) comprising (for example) porosity, V.sub.clay (clay content) and S.sub.w (water saturation) is used to calculate the equivalent geophysical properties. This is done using rock physics relationships that link the rock and fluid properties of the rock and fluid model to elastic properties or attributes (P-velocity and S-velocity, density, and impedance) (8) of the elastic model and also electric properties of the rock and fluid model (anisotropic resistivity) (9).

(35) Populating the rock and fluid model with rock and fluid properties may be done by using deterministic relationships, such as the soft sand or stiff sand correlations (for the elastic case) and Archie's equation or the Simandoux equation for the electric case.

(36) The correlation providing the relationships used should be calibrated to well log data. Alternatively, the rock physics models could be statistical relationships between rock and fluid properties and geophysical attributes, again based on well log data.

(37) Preferably the rock physics relationships should account for electric and elastic anisotropy in the sub-surface.

(38) Rock physics correlations according to another embodiment may also be facies dependent and therefore variable within the well, section or volume under study, according to the seismically defined facies that are observed.

(39) On the electrical side, the resistivity properties derived are then used to calculate the equivalent CSEM amplitude and phase response through a process of forward modelling of the 1D, 2D or 3D resistivity model that has been derived (10). Forward modelling algorithms to accomplish this task are well-known.

(40) In step (11) the elastic properties derived from the rock and/or fluid model are compared to the input data, the elastic model, and the misfit between the data and modelled attributes calculated.

(41) Similarly the CSEM response from the rock and/or fluid model is compared to the input CSEM data (12) and the misfit between the two calculated.

(42) If the misfit (either elastic or CSEM or a combination of the two) is higher than a user-defined level, the rock and/or fluid model is updated (13), and the process is repeated.

(43) The iterative update and comparison is repeated until the misfit between the input data and that modelled from rock and/or fluid model is less than a user defined value, or the maximum number of iterations is reached. The result is a final rock and/or fluid model (14).

(44) Several embodiments arise depending on the order of iteration over the seismic properties or the electromagnetic properties.

(45) In the second embodiment the inversion process identified within the dashed line as (6) in FIG. 3 is run using a cascade approach in which the elastic and electric processes are run sequentially rather than simultaneously within the general iterative loop. This is illustrated in FIG. 3 wherein now numbers start again from 1.

(46) As in the first embodiment, the input to the inversion process is elastic attributes in depth stored in an elastic model derived from the inversion of seismic data and, processed CSEM amplitude and phase data.

(47) In the first elastic iteration, starting from an initial rock and/or fluid model already populated with a first proposal of properties, elastic properties are calculated using elastic rock physics correlations or deterministic rock physics models as above (1), and said calculate properties compared with the elastic model. If the calculated misfit is high, the rock and/or fluid model is updated, and the comparison repeated, until a minimum misfit between the elastic properties calculated from the rock and/or fluid model and the elastic model obtained from the survey data.

(48) According to a specific example, the rock and fluid properties which are updated could be solely those that are best constrained by seismic data (for example porosity and clay content) or preferably comprises the complete set to be determined by the cascade inversion (for example porosity, clay content, V.sub.clay and water saturation, S.sub.w).

(49) Once a rock and fluid model that minimised the misfit to the elastic data has been found, this rock and fluid model is passed to the electric iterative component of the method, the electric iterative loop.

(50) Electric rock physics correlations or deterministic rock physics models (3) are used to calculate a resistivity model (which is anisotropic in this embodiment) (4), from which synthetic CSEM data are calculated using a forward modelling approach (5). These data are compared to the input CSEM amplitude and phase data. The rock and/or fluid model is then iteratively updated until a model that minimises the misfit between the CSEM modelled and measured data is found. As in the electric iteration, the update to rock and/or fluid properties could include all properties of interest, or be confined to those to which the CSEM data are primarily sensitive (for example water saturation).

(51) At the end of each elastic and electric iteration pair, the workflow is assessed for convergence (6).

(52) The convergence criterion could be based on the combined electric and elastic misfit falling below a user defined level, the change in the misfit between successive electric/elastic pass pairs falling below a user defined level, or a user defined maximum number of iterations being exceeded.

(53) FIGS. 4 (comprising FIGS. 4a and 4b) to 6 show an illustration according a specific case of the second embodiment.

(54) FIG. 4a shows a simple 1D model of porosity, clay content and water saturation (S.sub.w) wherein the value of each variable is represented in the abscissa axis and depth coordinate is represented in ordinate axis.

(55) From this synthetic elastic properties and CSEM data were generated. The elastic properties were generated using the stiff sand model and Gassman fluid substitution, a well-known rock physics approach (for instance disclosed in “The Hand book of Rock Physics” by Mavko et al.) linking elastic properties to the underlying rock and fluid properties. The resulting acoustic impedance and Poisson's ratio for the starting model are shown in FIG. 4b.

(56) In this example P-impedance (Ip) and Poisson's ratio (PR) were considered. Parameters in the stiff sand model were as follows: Pressure=1.0, critical porosity=0.4, coordination number=8, bulk modulus of brine, K.sub.brine=2.29, density of brine=1.0, density of hydrocarbon=0.065.

(57) Similarly the well-known Simandoux relationship (with parameters a=1.0, m=1.7, n=1.7, resistivity of clay, R.sub.clay=5.0 and resistivity of interstitial water, R.sub.w=0.18) was used to calculate the resistivity of the reservoir interval:

(58) $\frac{1}{R}=\frac{{\varphi }_e{S}_w^n}{{\mathrm{aR}}_w}+\frac{V_{\mathrm{clay}}{S}_w}{R_{\mathrm{clay}}}$

(59) where R is the effective resistivity of the medium, ϕ.sub.e is the effective porosity defined as ϕ.sub.total(1−V.sub.clay) being ϕ.sub.total the total porosity, V.sub.clay is the clay content, and S.sub.w is the water saturation.

(60) This model of reservoir properties was embedded in a simple 1D anisotropic background and from the resulting model the CSEM amplitude and phase response was calculated using a 1D forward modelling approach.

(61) A small amount of noise has been added to the data for the test. The CSEM amplitude and phase data, and the elastic properties calculated in this way form the input to the workflow proposed here.

(62) FIG. 5 shows the cascade approach to the joint inversion. In this example four iterations were undertaken: two electric iterations identified as step ii) in the invention and two elastic iterations identified as step i) in the invention.

(63) In each panel of FIG. 5 the solid black line is the true rock and fluid model (the solution of the problem that should be approximated after convergence of the method), and the dots represent the current rock and fluid model at each stage of the inversion cascade. Porosity and V.sub.clay are assumed constant in the reservoir interval in the starting rock and fluid model (1), and the reservoir is assumed water wet. Again, numbers just identify the sequence of steps.

(64) The first iteration is an elastic iteration (2), and for the example shown only porosity and V.sub.clay are updated in this step.

(65) A rock and fluid model is found which minimises the misfit between the modelled P-impedance and Poisson's ratio, and the input synthetic data calculated by simulating the response of the domain using a global search approach. Although the rock and fluid model is updated, it is still not close to the true rock and fluid model.

(66) The second iteration is an electric iteration (3), and in this example only S.sub.w (water saturation) is updated in this iteration.

(67) A rock and fluid model is found which minimises the difference between the modelled CSEM amplitude and phase data and the input synthetic amplitude and phase obtained by simulating the response of the domain when being excited with the dipole antenna, again using a global search approach. The updated S.sub.w is now less than 1 (water wet), indicating the presence of hydrocarbons in the reservoir.

(68) The third iteration is again an elastic iteration, and at this step, with the updated S.sub.w closer to reality, we find porosity and V.sub.clay, which are varied in this step, are now much closer to their true values.

(69) The fourth iteration is again an electric iteration (4), and this step again updates only the S.sub.w property, refining the value so that is now very close to the average value of the true S.sub.w in the reservoir.

(70) FIG. 6 shows the evolution of the misfit of the elastic properties and CSEM properties throughout the iterative process until convergence. This shows 6 iterations of the inversion described above, and it is clear that by iteration 4 the inversion has converged and there is little further change to the rock and fluid model. The final rock and fluid model accurately recovers porosity, V.sub.clay and the average water saturation in the reservoir interval.

(71) In the third embodiment the joint inversion is undertaken using an ensemble based method for the joint inversion of the electromagnetic data and seismic properties according to the invention to estimate rock and fluid properties, for example porosity, lithology, and saturation. The method according to this embodiment is illustrated in FIG. 7a which ends at FIG. 7b.

(72) Steps a) and b) are commonly carried out by any disclosed embodiment and therefore new embodiments emphasize the specific manner of generating the rock and/or fluid model and the steps iterating on it. In this specific case the rock and fluid model comprises a set of N.sub.e individual rock and fluid models m.sub.j, j=1 . . . N.sub.e (rather than just one) described by the statistics of the rock and fluid properties in the domain.

(73) As the rock and fluid model comprises a plurality of individual rock and fluid models each with different rock and fluid property attributes, each individual rock and fluid model m.sub.j is generated at the resolution of the acoustic model or the elastic model and falls within the range of possible rock and fluid properties consistent with the acoustic and/or elastic model. An initial set (the prior ensemble) of rock and fluid models are generated (in our case, models of porosity, lithology, and saturation) by populating the already generated acoustic or elastic model as in any of the previous embodiments providing individual rock and fluid models having different rock and fluid properties.

(74) The iterative method also computes the physical response (in our case, the elastic and electromagnetic data according to steps i) and ii) respectively), and all the individual rock and fluid models are updated using a Bayesian updating approach by minimizing the mismatch between the measured data and predicted response.

(75) The input data to the iterative loop is as in previous embodiments: the acoustic or elastic properties gathered in the acoustic or elastic model derived from seismic data, and processed CSEM amplitude and phase data obtained in the CSEM electromagnetic survey.

(76) The method operates as follows with reference to FIGS. 7a and 7b wherein numbers in circles identify again the order of the steps: 1) If the inversion is performed assuming a small error in the populated data, the risk is that we force the rock and fluid model to match the data by assuming unreasonable values (the solution is a local minimum but does not have a physical meaning). If the rock and fluid model is updated assuming a large error in the data, the risk is that the initial models are not perturbed enough (the solution is close to the initial model but the data match is poor). Therefore the updating process is repeated several times (number of updates or assimilations, N.sub.a, is typically between 4 and 8, and is chosen by the user in step (1)), starting with large errors, and reducing the error at each assimilation (5). Therefore the Bayesian updating is repeated N.sub.a times, as a stochastic optimization loop. The optimization is controlled by user defined parameters α.sub.i, which vary with assimilation, i, for i=1 . . . N.sub.a 2) An ensemble of initial realisations of the property or properties of interest is then generated (for example resistivity, porosity and/or saturation). The ensemble contains N.sub.e individual rock and fluid models (for example N.sub.e=100) (2) 3) For each of these N.sub.e individual rock and fluid models electric and elastic rock physics models or correlations are used as previously described to compute resistivity and elastic properties (for example impedance) (3) 4) The resistivity computed from each individual rock and fluid model allows an electromagnetic simulation to generate CSEM amplitude and phase data (4). 5) Based on the misfit between the observed and the calculated elastic properties and CSEM electromagnetic data, a Bayesian update is performed for each individual rock and fluid model comprising the ensemble (5) [identified in a general form as the rock and fluid model]. The Bayesian update requires the covariance matrix of the data and the cross covariance matrix of the rock and fluid model and data. These matrices cannot be analytically computed because of the non-linearity of the physical model. Hence said matrices are approximated using the ensemble of the individual rock and fluid models. The individual rock and fluid model update is performed as follows: a. For each individual rock and fluid model the vector of observed data is perturbed, as follows:
d.sub.uc=d.sub.obs+√{square root over (α.sub.i)}C.sub.D.sup.1/2z.sub.d where z.sub.d˜N(0,I.sub.N.sub.d) wherein d.sub.uc is the updated vector of observed data, d.sub.obs the observed property, that is, the elastic/acoustic properties determined from the seismic survey data and the electromagnetic data observed in the CSEM electromagnetic survey data, α.sub.i a number of perturbing updates being α.sub.i∈[1, N.sub.a] and, C.sub.D the covariance matrix, z.sub.d a sampled value of a Gaussian distribution. and then b. The updates of all the individual rock and fluid models m.sub.j comprising the ensemble of models, is performed using the following relationship all the:
m.sup.p:=m.sub.j.sup.f+C.sub.MD.sup.f(C.sub.D.sup.f+α.sub.iC.sub.d).sup.−1(d.sub.uc,j−d.sub.j.sup.f) wherein “:=” denotes the updating process as an update from former values, the superscript “p” denotes the updated value of the property, the superscript “f” denotes the previous model value of the property, C.sub.MD.sup.f is the cross covariance matrix between the prior vector of model parameters m.sup.f and the vector of predicted data d.sup.f, C.sub.DD.sup.f is the autocovariance matrix of predicted data, α.sub.i a number of updates being α.sub.i∈[1, N.sub.a] with the constraint Σ.sub.i=1.sup.N.sup.αα.sub.i=1 and, d.sub.j.sup.f is the predicted response from the m individual rock and fluid model, and d.sub.uc is as given above 6) Repeating steps 3-5 for each of the N.sub.a assimilation steps (6). 7) The final predicted rock and fluid model is the average of the updated individual rock and fluid models in the final ensemble. The uncertainty in the result is also obtained from the distribution of individual rock and fluid models in the final ensemble.

(77) The advantage of this approach compared to more conventional stochastic optimisation methods (for example Markov chain Monte Carlo (McMC) methods, which have been employed for CSEM inversion by e.g. Buland & Kolbjornsen (2012)) is its relative computational efficiency.

(78) The computational cost of McMC methods is generally relatively high and the method is not highly parallelizable, since the model solution at each iteration depends on the solution model at the previous iteration; therefore, the applicability of the McMC method to the joint inversion problem in 2D or 3D is unfeasible at the moment.

(79) In the method according to this embodiment, a set of individual rock and fluid models is run in each iteration; hence each iteration step can be parallelized for computational efficiency. A further advantage of this embodiment of the invention is that an estimate of model uncertainty is obtained.

(80) A simple example of the steps of the method is shown in FIGS. 8a, 8b and 9a and 9b. FIG. 8a shows a synthetic CSEM amplitude and phase dataset, generated by forward modelling from a well log based resistivity model (crosses) being said synthetic CSEM amplitude and phase dataset interpreted as data obtained from a CSEM survey allowing to compare the method according to this embodiment with said synthetic data. Solid circles show the response of the final rock and fluid model after reaching the convergence.

(81) Amplitude and phase at three frequencies (0.2, 0.5 and 1 Hz) are plotted against source-receiver separation of the CSEM survey.

(82) FIG. 8b shows synthetic elastic properties, in this case P-impedance and Poisson's ratio, in the reservoir interval of interest, generated from the same well as the CSEM electromagnetic data using an elastic rock physics model correlation (solid curve). In this case, synthetic elastic properties are also generated by a computer identifying said synthetic elastic properties with those elastic properties that would be obtained from the inversion of seismic survey data. This synthetic data allows the convergence of the method according to this embodiment to be validated. FIG. 8b shows the synthetic elastic properties using a solid line and the response of the final rock and fluid model when reaching the convergence criteria using crosses.

(83) The goal is to invert these two datasets using the method according to the third embodiment, to recover the underlying porosity, clay content and water saturation.

(84) In this example, the ensemble consisted of 1000 individual rock and fluid models (N.sub.e=1000), the number of assimilations or updates, N.sub.a, was 4, and the optimisation parameters, i.e. the elements of the matrix C.sub.D, were set to be equal to 0 outside the diagonal and 1/10 of the variance of the data on the diagonal.

(85) FIGS. 9a and 9b show the results of the method according to this third embodiment. FIG. 9a shows the mean of the final ensemble of the rock and fluid model, the one obtained after convergence of the method (grey curve), compared to the true model (black curve).

(86) The CSEM response of the final rock and fluid model is compared to the data (solid circles) in FIG. 8a. Similarly, the elastic response of the final rock and fluid model is shown in FIG. 8b (crosses).

(87) FIG. 9b shows the complete updated ensemble of 1000 individual rock and fluid models along with the true model (solid and thick black curve), with the spread in values indicative of the uncertainty in the estimation of porosity, clay content and water saturation (using a plurality of fine lines, one per individual rock and fluid model).

(88) A good agreement between the inverted porosity, clay content and saturation and the true values is achieved.

(89) The fourth embodiment of the invention uses stochastic sampling of reservoir properties derived from seismic data, followed by rejection sampling.

(90) Steps a), b) and c) are carried out providing the elastic model and then, a plurality of individual rock and fluid models with the resolution of the seismic model are generated populating the porosity, lithology and S.sub.w on the basis of the seismic model.

(91) The method according to this embodiment selects a sub-set among these individual rock and fluid models by matching their CSEM response to the measured data in the CSEM electromagnetic survey, retaining only those individual rock and fluid models with a misfit less than a user defined tolerance.

(92) This approach has two main advantages: firstly the resolution of the seismic is retained throughout the process, and secondly, this method provides an indication of the uncertainty in the result.

(93) The main steps according to this embodiment are illustrated in FIG. 10.

(94) In the first step of the method, seismic data, preferably pre-stack seismic data, are inverted to yield elastic impedance properties; that is, the elastic model (1). As in former embodiments, numbers only show the order of the method according to this example.

(95) Either a deterministic (such as the first and second embodiment) or stochastic (such as the third and fourth embodiment) inversion algorithm could be applied. If not already in depth, the resulting elastic properties must be converted from time to depth (2), to yield elastic properties of the elastic model in the depth domain (3).

(96) The elastic properties are then used to calculate rock and fluid properties (4) generating the rock and fluid model, for example porosity, clay content, water saturation (5).

(97) The resulting rock and fluid model may be 1D, 2D or 3D. Seismic facies may also be defined at this stage. Methods such as the multi-attribute rotation scheme (MARS) or statistical rock physics could be used for this purpose. In addition, probability density functions (PDFs) for the derived rock and fluid properties of the rock and fluid model must be calculated (6).

(98) If a stochastic seismic inversion process was applied in step 1 when generating the elastic model, then these can be derived from the inversion output. If not then PDFs can be derived from statistical rock physics approaches, or assumed based on prior knowledge of the likely distributions.

(99) The next stage is to sample the rock and fluid model using a stochastic approach (7). Note that the PDFs defined at different points in the elastic model may not be sampled independently, since the resulting models would be unreasonably rough and not representative of real geology. The realizations must therefore be spatially correlated. The realizations can be generated, for example, using the probability field simulation method of Srivastava (1992).
x.sub.i.sup.s(u)=x.sub.i.sup.mean(u)+σ.sub.i(u)∈.sub.i(u)

(100) where x is the property to be sampled (for example porosity, clay content or saturation), u(i,j,k) is the grid cell location (for instance in a 3D case) in the numerical model, x.sub.i.sup.s(u) is the simulated value at position u, x.sub.i.sup.mean is the mean value, σ.sub.i is the standard deviation of the property to be sampled, and ∈.sub.i is a spatially correlated error field.

(101) The trend of the mean can be obtained from the seismically derived inversion results, and the standard deviation either assumed, or calculated from the seismic inversion output (1) or the rock and fluid model (4).

(102) The number of individual rock and fluid models sampled, N.sub.e, is defined by the user, and could be for example 1000.

(103) For each of the N.sub.e individual rock and fluid models comprising lithology and fluid properties, N.sub.e individual rock and fluid models with electric rock physics properties are also generated in the rock and fluid model to calculate N.sub.e equivalent models of resistivity (8).

(104) Preferably these individual rock and fluid models with electric rock physics properties should account for electrical anisotropy in the earth so that the resulting resistivity models are also anisotropic.

(105) For each of the N.sub.e individual rock and fluid models of resistivity a forward modelling approach in 1D, 2D or 3D (depending on the dimensionality of the models) is used to calculate the equivalent CSEM amplitude and phase response (9) by simulation.

(106) In the final step the misfit between the observed and calculated responses are calculated (10). The misfit could be defined as (for example) the .sup.2-norm between calculated and observed values.

(107) Any model with a misfit greater than the user-defined tolerance is rejected, leaving a population of N.sup.accepted final models. Alternatively the user could define N.sup.accepted, and accept the N.sup.accepted models with lowest misfit (11).

(108) FIGS. 11a, 11b and 12a and 12b show a simple synthetic example of the proposed method.

(109) FIG. 11a shows the porosity, saturation and clay content within a reservoir (thick black solid lines). In the proposed method these would be estimated using seismic data.

(110) In this example porosity and clay content are assumed known. Therefore 100 realizations of the water saturation have been generated (FIG. 11a, fine lines).

(111) The simulations are generated using a vertical correlation range (50 meters) and assuming an average saturation of 0.5 in the entire interval.

(112) The Simandoux model, which links electrical resistivity to the underlying porosity, clay content and water saturation was then used to generate 100 equivalent resistivity models, and the CSEM amplitude and phase response of each of the models was calculated using a 1D forward modelling algorithm. The resulting responses are shown using fine lines in FIG. 12a. In this figure the CSEM amplitude and phase of the true model (the observed data collected in the CSEM electromagnetic survey) is shown by the solid circles, for signal frequencies of 0.2, 0.5 and 1 Hz.

(113) FIG. 12b shows the misfit between the observed CSEM response and the response calculated from the 100 model realisations. In this example the 10 realisations with the lowest misfit were accepted: these are plotted in FIG. 11b. Although there is still some variability in the responses, the average saturation in the reservoir interval (which is well constrained by the CSEM method) is recovered accurately.