• Bayesian inversion
• carbonate reservoir
• elastic inversion
• rock physics model
• time-lapse

When simulating a production scenario of a reservoir, it is fundamental to have a rock physics model, that allows to predict the variation of subsurface elastic properties, during the different production stages. The time lapse elastic inversion technique, aims to retrieve from seismic data, the variation in time of the elastic properties describing the subsurface. In this context, the presented project, aims to apply a Bayesian approach in order to solve a time lapse inversion problem. Furthermore, an analysis of the solution uncertainties and stability are also presented, for different initial conditions such as different noise level in the seismic data set, or the case in which the frequency of the ricker wavelet used in the forward operator during the inversion process doesn’t match the source wavelet frequency.
The thesis belongs to the framework of a reservoir characterization process, and can be summarised in two stages: the first part, included the implementation of a rock physics model in a non clastic reservoir, in order to estimate the behaviour of the subsurface elastic properties for various reservoir fluid saturation conditions. The second part, consisted in exploiting such predicted properties, in order to compute the observed data for the selected reservoir saturation conditions, and hence applying an ensemble based Bayesian inversion in order to evaluate the subsurface elastic properties and their associated uncertainties at different time.
More in depth, during the first phase of the project four well logs dataset belonging to the same reservoir have been provided, as well as information relative to the fluid in place, and reservoir conditions. The geological setting, is represented by an isolated carbonate platform. Therefore, being a carbonate environment, characterised by a strong heterogeneity, it was crucial to set up a rock physics model that allowed to differentiate the presence of different pore types and map their distribution in order to obtain more reliable results when performing the fluid substitution workflow using Gassmann equation. Such task has been successfully accomplished by combining the Xiu Payne method for RPM and the Eshelby-Walsh theory: through the use of the relationship between P sonic and effective porosity log curves, the model here presented, allows to discriminate between three different pore types (interconnected refence pores, stiff pores and microcracks) and assigns a specific pore aspect ratio to each of them. Furthermore, a modified Gassmann equation, (Elshby-Walsh theory) has been applied, which includes a rock-matrix term m, function of the pore aspect ratio distribution, hence enabling a more realistic modelling of the bulk rock elastic properties for different fluid saturation conditions. It must be remarked here, that in reservoir characterisation, the use of a reliable RPM is fundamental, in order to build a reliable static model of the subsurface otherwise, there may be strong misleading estimations when using such static model as starting point for the dynamic reservoir simulations. The fluid replacement workflow allowed to simulate five different reservoir saturation scenarios, going from the in situ conditions, 100% gas saturation, to the fully brine saturated reservoir, with a constant increment of 25% in brine saturation. Batzle Wang equations have been used to model gas and brine elastic properties, and finally, Reuss-Voigt-Hill average has been applied when estimating the density and bulk modulus of the two phases fluid. Having obtained the values for P velocity, S velocity and Density at each simulated scenario, the corresponding 3D reservoir models of the subsurface have been created, using Jason software. Such models, have been exported as seg-y files and used as starting point for the second part of the thesis.
In the second phase of the project, due to computational time, three scenarios have been selected among the five computed in the first phase, corresponding to 100% gas saturation, an intermediate scenario with mixed fluid saturation (50% gas and 50% brine), and finally 100% brine saturation. Hence, the corresponding seg-y files have been imported in Matlab, and used to compute the synthetic seismogram representing the observed data of the inverse problem. Through the use of a convolutional forward model, based on the Zoeppritz equations, the VP, VS, and density models have been used in order to compute the pre-stack seismic response. Finally, an ensemble based Bayesian inversion has been implemented. As suggested by its name, an ensemble of prior models has been generated and iteratively updated according to the perturbation of the observed data and the Kalman filter, which, in turn, is computed as function of the covariances of the observed and predicted data and the cross covariance between the prior model and the corresponding seismic response. The inversion strategy here proposed consists in the implementation of the ensemble based Bayesian inversion, independently for each saturation condition. This method has been tested not only in ideal situation, (low noise level in seismic data and a perfect reconstruction of the wavelet in the Forward model), but also for strong noise conditions and in situations where the frequency of the wavelet implemented in the forward model operator didn’t match the frequency of the wavelet used to construct the synthetic seismograms (observed data).
The final step of the second stage, consists in the results analysis, where the misfit between the elastic properties obtained from the inversion and the elastic properties modelled in the first stage of the project, is evaluated in order to assess the stability of the inversion algorithms. The estimated parameters are in accordance with the modelled once, thus, demonstrating the reliability of the implementation of the Ensemble Based Bayesian inversion. A dramatic increase of the misfit has been noticed when the observed data are contaminated with higher level of noise. Furthermore, it has been noticed that the elastic properties of the subsurface aren’t properly recovered, showing unreliability of the inversion process, in the case when the frequency of the Ricker wavelet used in the forward operator is lower than the frequency of Ricker wavelet used to create the observed data. On the contrary, results are still acceptable if a Ricker wavelet with a frequency higher than the one used for the creation of the synthetic data is implemented in the forward operator.
A further development of the project could be represented by the implementation of the Discrete Cosine transform (DCT) or Singular Value Decomposition (SVD), in order to compress the size of the model and data space. Such procedures would allow to speed up the computational time of the inversion, and hence decrease the computational cost. An additional strategy, could be represented by the implementation of a joint time lapse inversion.