Abstract

Offshore gravimetric monitoring has been introduced as a complement to seismic monitoring of fields with moving fluids. The Sleipner field in the North Sea is a fully operational carbon capture and storage facility, where CO2 is injected for storage. Gravimetric measurements are one of the geophysical monitoring methods applied, and the data have been used to estimate the in situ density and dissolution of the CO2. We defined a Bayesian inversion of gravimetric data, and we used this to analyze gravimetric data at Sleipner field. In our approach, we included spatial uncertainty in the model and performed a Bayesian analysis of the in situ CO2 density and dissolution. We also analyzed the impact of mass changes due to gas production from the Ty Formation. Our estimates were comparable with published results.

Introduction

The Sleipner field in the North Sea is a fully operational carbon capture and storage facility. The CO2 is injected in a saline aquifer — the Utsira Sand, and the migration of the gas is monitored by time-lapse seismic surveys, time-lapse seabed gravimetric measurements, electromagnetic surveys, and seabed mappings. The key aims of the monitoring are to ensure safe storage of CO2 within the storage reservoir, to keep track of the distribution and migration of the gas through the reservoir and possibly into adjacent strata, and for early detection of leakage toward the seabed (Chadwick et al., 2006). Interpretation of time-lapse seismic surveys gives a good image of the distribution and migration of the CO2 plume, including the depth, whereas gravimetric measurements provide complementary constraints on the in situ density of the injected CO2.

For the Sleipner field, the first gravimetric measurement in August 2002 serves as the baseline for time-lapse gravimetric surveys. As CO2 is injected into the storage reservoir, it fills the pore space, displaces the brine, and reduces the total density of the rock. The density decrease has an effect on the local strength of gravity. Monitoring these changes in the gravity strength can give constraints on the in situ density of CO2. For the Sleipner field, there is a significant uncertainty in the storage reservoir condition regarding temperature and pressure for all stages in the process, that is, before, during, and after injection. At the predicted temperature and pressure conditions, the CO2 is close to the critical point of phase transition. This means that small changes in reservoir temperature may cause large changes in the CO2 density. The spatial variability, which has not been incorporated in the analysis previously, might give important contributions to the uncertainty.

Previous work by Nooner et al. (2007) uses time-lapse gravity measurements together with forward gravity modeling based on seismic data and reservoir simulation models to constrain the in situ density of CO2 within the Utsira Formation. They provide a best-fit average in situ CO2 density of approximately 530±65kg/m3 (95% confidence interval). Later work by Alnes et al. (2008) has reprocessed the original gravity data, resulting in a new best-fit average density of CO2 of 760kg/m3, with an upper bound of 770kg/m3 and a lower bound of 640kg/m3 at 95% confidence. Both of these works base their analysis on time-lapse data from 2002 to 2005. Alnes et al. (2011) include gravimetric measurements from 2009 and provide an updated best-fit average density of 720±80kg/m3, assuming there is no dissolution of CO2 in the formation water. They also combine the gravimetric data with improved temperature measurements to provide an upper bound on the dissolution of CO2 in formation brine.

Bayesian inversion of geophysical data contributes to the analysis of geophysical data by mapping the uncertainty in data into uncertainty about earth parameters (Tarantola and Valette, 1982; Duijndam, 1988a, 1988b; Scales and Tenorio, 2001; Malinverno and Briggs, 2004; Bosch et al., 2007). Several types of geophysical data with different sets of models have been studied within this context (Gouveia and Scales, 1998; Buland et al., 2003; Gunning and Glinsky, 2004; Hoversten et al., 2006; Buland et al., 2011). Rock-physics models combined with geophysical data can be used to obtain a quantitative prediction of properties in the subsurface reservoirs, e.g., Mukerji et al. (2001). It is common to use a two-step approach by first performing a Bayesian inversion of the geophysical data and next interpret the model in light of the uncertainties in the inversion (Doyen, 2007; Buland et al., 2008). Hauge and Kolbjørnsen (2014) present the Bayesian inversion of gravimetric data and show quantitatively the information content held by the data in terms of which features that are resolved.

In this paper, we use the two-step approach to geophysical inversion; that is, we define a Bayesian inversion of gravimetric data and show how this can be used for analyzing properties derived from gravimetric data. In particular, we do a Bayesian analysis of gravimetric data from the Sleipner project. We show how a geostatistical framework makes it possible to account for the spatial uncertainty and include multiple sources of information. We also discuss how the Bayesian approach can be used to analyze the contributions due to gas production at the Ty Formation. Our estimates for in situ CO2 density and dissolution at the Sleipner injection site are in good agreement with estimates provided in the existing literature.

Method

We present here the Bayesian inversion of gravimetric data using a geostatistical methodology, and we show how to include the contributions from an additional source of mass changes.

Bayesian inversion of gravimetric data

In the Bayesian approach, it is required that a priori knowledge about the field under study is quantified. In the field of reservoir characterization, there is substantial knowledge about modeling of rock-physics properties at the geomodel scale (Dubrule, 2003; Eidsvik et al., 2004; Doyen, 2007). This is, therefore, a natural scale to build the a priori information into the model. This model is, however, too fine to make efficient computations; thus, we perform an intermediate step of upscaling to provide the model with uncertainty on a computational/coarse scale grid. To assess the error in the upscaling, we generated 100 models from the a priori distribution and computed the gravimetric response using the fine-scale and the upscaled approaches. The increase in the uncertainty added in the error standard deviation when including this component is less than 1%. The coarse-scale grid is thus more than sufficient for resolving the information content in the gravimetric data. By introducing the geomodel as a reference, we are able to include rock physical and geostatistical knowledge such that the variability at the coarse-scale model is physically based.

Gravimetric data are a truly 4D data type in the sense that the meaningful content is the change in gravimetric response between two vintages. We thus consider the density change Δρ(x) as the relevant earth model parameter, with x as the spatial position. We model the parameter as a Gaussian random field m(x)=Δρ(x):  
m(x)=N(μm,Σm).
(1)
Here, μm is the expectation of the model parameter and Σm defines the covariance structure. In geostatistical inversion, this covariance imposes spatial dependencies in the model. We will discuss how this model is built for the Sleipner CO2 injection below.
For gravimetric inversion, Newton’s law of universal gravitation provides a linear relation between the earth parameters and the data. We therefore have a linear relation between the change in gravimetric data d and the change in density m:  
d=Gm+ϵ.
(2)
The matrix G has entries that describe the vertical component of the gravimetric contribution of each position in the inversion region to all measurement positions:  
Gi,j=γ|dx|xixj2cosθz.
(3)
Here, γ is the gravitational constant; the grid cells have volume |dx|; and the spatial positions xi are the positions in the inversion region, looped through by index i, whereas index j loops through the measurement positions at the sea bed. For the vertical component, we multiply with cosθz, where θz is the angle between the full attraction and the vertical axis. In other words, the elements of matrix G are inversely proportional to the square distance between measurement positions and all positions in the modeled reservoir, thus inherently a smooth matrix. We consider changes of gravimetric response in time due to density changes in the reservoir, and we assume that outside contributions are accounted for in the processing of the data (Nooner et al., 2007; Alnes et al., 2008; Zumberge et al., 2008). After the processing, an error in the observations still remains and is denoted by ϵ in equation 2. The error is modeled to have two terms:  
ϵ=ϵw+1ϵc.
(4)
A general white-noise error term ϵw accounts for independent measurement errors given as  
ϵw=N(0,σw2I),
(5)
with I being the identity matrix. Finally, a scalar error term ϵc, common for all data, is added to the linear system to include the possibility to correct for shifts in the level of data.
The only changes we solve for are in the inversion region we define. We point out that multiple change regions may be included in the inversion as discussed below. The a posteriori distribution for the model parameter given the gravimetric data, m(x)|d=N(μm|d,Σm|d), is given by standard results as  
μm|d=μm+ΣmGT(GΣmGT+Σϵ)1(dGμm),Σm|d=ΣmΣmGT(GΣmGT+Σϵ)1GΣm.
(6)
The white noise and the common error term are incorporated in Σϵ=σw2I+σc21T1, with the vector 1 having ones in all entries and having length equal to the number of measurement positions.
A central part in computing the a posteriori distribution is the expression GΣmGT that involves multiplication of a matrix of dimensions equaling the number of grid cells in the inversion region. For a geomodel, the number of grid cells is potentially very large, following that the multiplication is potentially computationally very expensive. To handle the large size, the conditional inverse problem is upscaled to a coarser reservoir with fewer grid blocks (the computational model). Here, we introduce the term block for the collection of fine-grid cells that make up a grid cell in the coarser computation grid. Only quantities related to the inversion grid are upscaled, that is, the a priori distribution described by μm and Σm and the gravity matrix G. Upscaling the latter matrix follows straightforward upscaling based on coarsening the grid quantities. As described above, the G matrix is a smooth, low-frequency matrix, and for all upscaling relevant for this paper, the error due to the upscaling was several magnitudes less than the observation errors. Thus, the upscaling of this matrix does not introduce any significant loss of accuracy. The stochastic variable m is upscaled by linear block averaging given by the following formulas:  
mup(A)=1NAxiAm(xi)
(7)
and  
Cov(mup(A),mup(B))=1NANBxiAxjBCov(m(xi),m(xj)),
(8)
where A and B denote grid blocks. The computations in equation 8 can be done efficiently using the standard convolution trick through the Fourier domain. We point out that the analyses of the gravimetric inversion in the following sections are performed on the upscaled computational model. Under sufficient regularity assumptions for the sampling design, it is possible to bring the inversion result back to the geomodel. This is, however, not relevant for the discussions below, and we leave it out.

Including additional reservoirs in inversion model

We will perform the Bayesian analysis of data from the Sleipner injection project. These data are influenced by gas production from the Ty Formation, which is located laterally close to the Utsira injection site. To incorporate this in our analysis, we need to introduce an additional reservoir with density changes.

First, we repeat the mathematical formulation of the Bayesian inversion setup, now with an additional parameter modeling the additional mass change. The contribution from an additional reservoir to the data d is given through an additional term in the linear relationship:  
d=G1m1+G2m2+ϵ.
(9)
Now, G and m defined previously in the section on the methodology are redefined for our model to have subscript 1, that is, G1 and m1. For the additional reservoir, we use G2 and m2 for the gravitational model matrix and the earth parameter, respectively. Similarly, the additional reservoir has a covariance matrix Σm2 describing the spatial correlations for this reservoir. The error term ϵ is the sum of measurement error and the common scalar error term to adjust the level of data, as defined previously. It would be possible to solve for both density changes simultaneously, but our primary interest is related to the first reservoir. The new a posteriori distribution for the reservoir of interest is now given as  
μm1|d=μm1+Σm1G1T(G1Σm1G1T+G2Σm2G2T+Σϵ)1(dG1μm1G2μm2),Σm1|d=Σm1Σm1G1T(G1Σm1G1T+G2Σm2G2T+Σϵ)1G1Σm1.
(10)
By comparison to the expressions in equation 6, the a posteriori distribution is altered by having an additional term G2Σm2G2T parallel to the error matrix expression, as well as G2μm2 parallel to the data mean. Thus, in our analysis, the contribution from the additional reservoir gives a bias correction and an additional error term consisting of coherent noise.

Bayesian analysis of inversion results

In our analysis, the main interest in the inversion is total mass change. This quantity can be related to the in situ CO2 density and dissolution (Alnes et al., 2011). The total mass change in the initial reservoir can be deduced from the density change:  
ΔM=DΔρ(x)dx.
(11)
For our model parameters, this is given by the relation:  
ΔM=|dx|×1Tm.
(12)
This gives the a posteriori distribution of total mass change as 
p(ΔM|d)=N(|dx|×1Tμm1|d,|dx|2×1TΣm1|d1).
(13)
Reformulating slightly the notation in Alnes et al. (2011), we find that ΔM, the mass change; ρCO2, the average density of CO2; ρB, the average density of displaced brine; MCO2, the mass of injected CO2; and α, the dissolution have the following relation:  
ΔM=MCO2(1α×k)(1ρBρCO2)+MCO2×α×k,
(14)
where k is a constant specific for the injection regime. The first term compensates for displacement of brine by free CO2, and the latter term is the contribution of the dissolved CO2 to an increased brine density. Thus, we assume that the volume of brine does not change when absorbing CO2.

Bayesian analysis of the Sleipner gravimetric data

The gravimetric data analyzed here are from the Sleipner injection site. We discuss the 2002–2005 and 2002–2009 changes. The data have been processed to remove the effects of the underlying Ty Formation, which is a gas-producing field deep below the Utsira formation and located laterally very close to the injection site (Alnes et al., 2008, 2011).

Defining statistical models

We consider analysis with and without accounting for the uncertainty in the Ty Formation. In addition, we do the analysis with an open a priori distribution and an a priori distribution constrained by pushdown from seismic data. The pushdown is the change in time interpretation of a seismic reflector below the reservoir when comparing two surveys. Substitution of CO2 into the rock reduces the seismic velocity and causes the reflector to appear deeper on the monitor survey. An overview of the inversion cases is given in Table 1.

We carry out gravimetric inversion with the Sleipner data on a storage reservoir extending 2500 m in the east direction and 5500 m in the north direction and with thickness of 200 m. The top of the reservoir is placed 820 m below mean sea level. The geomodel grid is resolved into 100×220×100 cells with cells of size 25×25×2m3. The computational grid that we perform the inversion on is resolved into 10×22×4 blocks with size 250×250×50m3. The lateral geometry of the data and inversion region is given in Figure 1. When building the a priori model, information from alternative sources might be used to provide a rationale for the assumptions.

According to Alnes et al. (2011), temperature considerations suggest that the average CO2 density in the Utsira region is 0.675±0.02g/ccm. Using the center of this interval and if we in addition neglect dissolution and use the injected mass in Table 2, we can provide a rough estimate for the mass change using equation 14. In the basic a priori distribution, we divide this equally in all grid cells in the inversion region. For the a priori distribution constrained by seismic, we divide the mass evenly along a vertical profile. Because the pushdown gives a strong indication of the presence of CO2, we reduce the uncertainty in density outside the pushdown region. Moreover, because the magnitude of the pushdown is an indication of the amount of CO2 in a profile of the CO2, we distribute the mass changes proportional to the pushdown change, i.e., the pushdown from 2001 to 2006 and 2001 to 2008, respectively. Note that because we allow full variability around this mean, the a posteriori density change will not be proportional to the pushdown. Figure 2 shows a smoothed version of the pushdown for the years 2001, 2006, and 2008, respectively.

The mean density change is only one aspect of the description of the density changes. The uncertainty in the geomodel is described by a covariance function. In the basic a priori distribution, we use an exponential covariance function, which has the following form:  
Cov(m(x),m(x+Δx))=σm2×exp(3(ΔxRx)2+(ΔyRy)2+(ΔzRz)2).
(15)
The parameters into the expression are the standard deviation of the local density change σm and the ranges (Rx,Ry,Rz), which describe the extent of typical features in the density change. To assess the variability at the geoscale, we consider Figure 3, which illustrates the full uncertainty span in brine and CO2 density. The difference in the density will typically be less than 0.550g/ccm. In addition, local changes in saturation and porosity add to the spatial variability to how large the changes might be. According to Arts et al. (2004), the porosity is measured as large as 42% and the saturation might become as large as 90%. Together, this leads to local density changes of the order of 0.2g/ccm. The uncertainty in the geomodel should cover this value range, which suggests the standard deviation σm=0.1g/ccm. From the pushdown, we see that the size of typical features in the difference is of the order 500 m in both lateral directions. In the vertical direction, we set the range to be 10 m, which is about the thickness of the CO2 layers. In the a priori model constrained by the seismic, we set the variance outside the region with pushdown changes to be small, i.e., σ=0.01, because there is not likely to be any saturation changes in these regions. In addition, due to the layered structure of the CO2 indicated on seismic amplitudes, we assume that about half of the layers only have small saturation changes. This will further reduce the variance in the upscaled model with a factor of 2. The mean and standard deviation of the vertical sum in the upscaled a priori distributions are displayed in Figure 4.

The Ty Formation is located approximately 2300 m below mean sea level. The lateral region with the disturbance is approximately described by a square stretching from the southwestern point with Universal Transverse Mercator (UTM) coordinates (434, 6466) km and 3.5 km to the east and 8.5 km to the north; see Nooner et al. (2007). The thickness of the zone is very small compared to the resolution in gravimetric data; thus, we model only the spatial distribution of lateral changes. The ranges are set to be 500 m in both lateral directions. The gravimetric data analyzed in this section were processed to remove the effect of the Ty Formation (Alnes et al., 2008, 2011). Therefore, we consider the bias in equation 10 to be removed, but the uncertainty remains; thus, only the covariance is required. Moreover, Alnes et al. (2008) estimate that the total mass has an uncertainty in the order of 5 MT. To give room for spatial deviations from the mean scale, the uncertainty is set such that the total mass change has a standard deviation of 5 MT. In a cell of the upscaled grid of 500×500m, the uncertainty is 0.26 MT, which corresponds to 1050kg/m3.

Gravimetric inversion

The measurement positions relative to the lateral outline of the reservoir are shown in Figure 1. We will consider two sets of measurements: The first set of data is the change in gravimetric response from 2002 and 2005, and the second set is from 2002 to 2009. The standard deviation for the measurement errors (given as part of the data) σw is approximately 3μGal, and the standard deviation for the scalar level shift is set to σc=1μGal. Figure 5a and 5b shows the processed data (Alnes et al., 2008, 2011) together with the response from the a priori mean as defined above.

The results of the Utsira inversions are displayed in Figures 6 and 7 for both data sets. The corresponding data fits are displayed in Figure 8a and 8b.

It is interesting to notice that both inversions (see Figures 6 and 7) suggest positive changes slightly off target, but compared to the uncertainties in the standard deviations, these are not significant. For the basic a priori distribution, the northernmost region has been very little influenced by the data and the uncertainty remains almost at the size of the a priori distribution. There is very little visual influence when including the Ty Formation. Figure 9 shows the effect on data caused by including the Ty Formation. Inspecting the inversion of Ty in Figure 10, we find weak indications of a positive contribution in the northern part. The uncertainties in the inversion are, however, large, and the spatial resolution is very poor. This reflects the lacking ability to constrain the spatial model.

Bayesian analysis of in situ CO2 density and dissolution

The target is to analyze the CO2 density and dissolution constrained by gravimetric data. Equation 14 relates this to the mass change. Table 3 shows the estimated mass changes and corresponding uncertainties. Note that all a posteriori standard deviations are considerably smaller than the a priori ones; thus, we might conclude that the gravimetric data carry much information about the total mass change. For the unconstrained a priori distributions, the a posteriori uncertainty is still much larger than the estimated effect. Thus, for these cases, we are not able to resolve the mass changes to a significant degree. In general, we see that inclusion of the Ty Formation results in a larger mass reduction; also, the uncertainty in the estimate is increased.

Alnes et al. (2011) report that k=10.0yr for the 2002–2009 survey. Assuming a constant rate of injection between 2005 and 2009, we find that k=7yr for the 2002–2005 survey. Assuming that the density is in the range 0.40.8g/ccm and that the dissolution is in the range 0.0%–10.0%, we can compute the likelihood of the combinations. This is shown in Figure 11 for the seismic constrained cases. We find that 2009 data are better resolved and have the mode of the mean density at higher values than the 2005 data. This reflects the larger uncertainty in the 2005 data. It is also reasonable that the 2005 data could have a lower average density if a larger fraction of the injected CO2 has cooled down in 2009. Note that the data do not separate the two effects, that is, density and dissolution. If the focus is on the dissolution and we assume that the average density of CO2 from both survey times is independent and has a normal distribution centered in 0.675 and a standard deviation of 0.01, according to computations from Alnes et al. (2011), we can compute the a posteriori distribution for the dissolution by summing out the densities. These a posteriori distributions are shown in Figure 12. The mean value is 0.67% if we neglect the effect of the Ty Formation. Including the Ty Formation in the analysis, the mean increases to 0.72%. In either case, the probability of exceeding 1.76% dissolution is less than 5%.

Conclusion

We have shown how it is possible to account for spatial variability in gravimetric inversion. For the Sleipner data, the lateral resolution is poor, but when we constrain our model with seismic information, we are able to make quantitative predictions about the CO2 density and dissolution. Using the restriction of the spatial average of the density provided by Alnes et al. (2011) and analyzing the 2009 data only, we find an upper limit (1.83%) for dissolution that is rounded down to the number provided by Alnes et al. (2011). Including the data from 2005 as well reduces the upper limit, but the value is still rounded up to the number provided by Alnes et al. (2011). The upper bound we obtain for dissolution is 1.76% per year, whereas the mean value is approximately 0.7% per year. The bound is sensitive to the assumptions related to the average CO2 density.

Acknowledgments

We thank the Research Council of Norway, Statoil, and ExxonMobil for their financial support for the project “Monitoring geological CO2 storage: Quantitative CO2 prediction with uncertainty from physical modelling and multiple time-lapse data types” under which part of this work has been done. We also thank the Sleipner license partners, Statoil, ExxonMobil, and Total for permission to publish Sleipner field gravimetric data, and finally, we thank A.-K. Furre, H. Alnes, and O. Eiken for useful discussions on this topic.

Biographies and photographs of the authors are not available.

Freely available online through the SEG open-access option.