Four Dimensional Gravity Forward Model in a Deep Reservoir

In this work, we calculate the gravity signature of small density changes in a real-case deep reservoir. Based on the 3D forward density modeling of constrained geometries and using parameters of the involved rocks and fluids, we compute the differential gravity signature before and after the production period of the Volve oil field in the North Sea. Causative sources of the retrieved residuals are spatially correlated with positions of the most productive wells, locating areas of maximum density change. Results show that the 4D gravity forward model is capable of resolving residual gravity signatures also for deep and small density changes. In particular, we locate ∼−13 μGal gravity minima over the 2750 m deep reservoir; this minima was caused by the −53 kg m–3 density change related to production and injection activity.


INTRODUCTION
Using gravity data for sub-surface modeling of geological structures may represent a useful tool to address open questions about geological or geophysical processes at the crustal or local scale. Application of 2D or 3D forward or inverse techniques depends on real gravity data availability over the target area (e.g., Mancinelli et al., 2015;Dressel et al., 2018;Fedi et al., 2018;Mancinelli et al., 2019Mancinelli et al., , 2020Sobh et al., 2019) but often the smaller gravity signatures are concealed and hard to elaborate through filtering. In the last two decades, the increase in gravimeters' precision allowed the development of 4D gravity models monitoring fluid-related density changes at increasing resolution (e.g., Hare et al., 1999;Eiken et al., 2000Eiken et al., , 2008Ferguson et al., 2007Ferguson et al., , 2008Vasilevskiy and Dashevsky, 2007;Alnes et al., 2008;Davis et al., 2008;Gasperikova and Hoversten, 2008;Stenvold et al., 2008;Jacob et al., 2010;Krahenbuhl et al., 2011;Krahenbuhl and Li, 2012;Wilson et al., 2012;Reitz et al., 2015;Braun, 2016, 2017). It should be noted that in all these case histories, the top of the gravity source was always above 2500 m depth. The depth of the source is the most critical parameter affecting the resolvability of a gravity anomaly (e.g., Blakely, 1996). In the case of CO 2 storage, the depth of 2500 m is considered a threshold between shallower (thus resolvable) and deeper (very difficult to resolve) field applications (Cooper et al., 2009).
The evolution of gravity modeling techniques has closely followed the technological improvements that allowed significant increase in precision and accuracy of gravimeters over the last decades. A compelling revision is provided by Van Camp et al. (2017). The goal of lowering the threshold for stable measurements, together with the increased ability of identifying noise sources, has allowed to reach precisions below 1 nm s −2 -i.e., 0.1 µGal, also for transportable absolute gravimeters (Ménoret et al., 2018). Seafloor gravimeters can achieve measurements with time-lapse accuracy better than 5 µGal (Eiken et al., 2000;Sasagawa et al., 2003;Stenvold et al., 2008). However, the monitoring of fluid-induced density changes at production or storage sites through gravity 4D models has not evolved in a similar way. In fact, very few production and storage reservoirs have been investigated using 4D gravity methods. This is likely due to problems related to three main sources: (i) size of the production/injection process in terms of depth, thickness, volumes, and densities of the involved deposits and fluids; (ii) high-quality gravity data availability due to costs for gravimeters and data acquisition and knowledge about noise sources; and (iii) best familiarity with the 4D seismic method which is preferred to any other approach. If the first two are related to the nonuniqueness of the modeling solution, the third is related to the results offered by modern 3D and 4D seismic surveys, even if the method has also limitations in locating fluid content, porosity, and saturation changes (Devriese and Oldenburg, 2014).
In this work, we calculate the gravity response of small density changes in a small and deep reservoir with a size/depth ratio slightly lower than 1. This scenario represents the worst case for the application of gravity methods to locate density changes related to fluid production or injection. Through 3D forward modeling of the reservoir geometries as obtained from seismic and borehole data, we model production-related density changes to retrieve the gravity response before and after production. Modeling is constrained by volumes and densities of produced and injected fluids at reservoir conditions.

The Volve Field Database
The Volve field is located in the Central North Sea (Figure 1A), and it was discovered in 1993 by the well 15/9-19 SR. It is located ∼5 km north of the Sleipner East field (e.g., Alnes et al., 2008Alnes et al., , 2011 and production started in February 2008 (Production data report in the Volve field dataset). Oil and gas production ended in September 2016 after a total production of 1.5 × 10 9 Sm 3 including oil, gas, and formation water. Among the several wells drilled in the field, the majority of the production (∼98%) was achieved by wells 15/9-11, 15/9-12, and 15/9-14 (Figure 1). The Volve reservoir is located in the Hugin Jurassic sandstones between 2750 and 3120 m total vertical depth (TVD) below sea level (Discovery report in the Volve field dataset). Hugin sands were deposited in a nearshore marine setting, show high total organic carbon (up to 80%), and are gas prone with the capability to generate also light liquid hydrocarbons (e.g., Isaksen et al., 2002). According to well data in the Volve field, the thickness of the Hugin formation can range between a few tens of meters up to ∼150 m where all the 18 levels of the formation have been drilled. Production at Volve was sustained by formation water injection of a total volume of 3 × 10 7 Sm 3 from the wells 15/9-F-4 and 15/9-F-5 that were active between April 2008 and September 2016 (Production data report in the Volve field dataset).
After shut down and removal of the production equipment in summer 2017, the entire dataset regarding seismic volumes, well logs, and petrophysical and geological characterization of the field was disclosed by Statoil (now Equinor) in June 2018 1 .
The Volve dataset encompasses over 5 Tbytes of data counting over 11000 files including both raw data and interpretations of the formation tops and faults. It includes also the 3D models used for production simulations as resulting from the seismic data interpretation and well logs. Locations of production and injection wells, daily and monthly reports on production, and injection activities are available as well. Finally, geochemical analyses of formation fluids at reservoir conditions of pressure and temperature are provided together with detailed characterization of the involved formations. The availability for the scientific community of such a detailed and complete database regarding a reservoir system and its production history represents a unique opportunity and allows the community for unprecedented attempts in testing modeling techniques.

Modeling Procedure
In order to provide an input geometry for the forward gravity model, we first created a 3D model of the reservoir using Hugin top and bottom surfaces as provided in the Volve dataset. The surfaces of the main faults in the area were depth-converted as they are provided as two-way travel-times (TWT) in the dataset only. This step was carried out using TWT-depth data provided in well logs and reports on well activities included in the dataset. Finally, the 3D model of the reservoir was obtained by intersection of the top and bottom of the Hugin formation with the bounding faults ( Figure 1B). We chose to not include other structures in the model to keep the geometry as simple as possible; moreover, we assume that all bounding faults are sealing -i.e., all density changes are contained within the reservoir volume.
Once the geometry of the reservoir was defined, we discretized the ∼4.3 × 10 9 m 3 Hugin volume in 275642 cubic cells with size 25 m × 25 m × 25 m and included it in a larger 3D mesh with same cell size. In this geometric model, the only surface included together with the Hugin volume is the sea bottom representing the major density contrast between sea water and the underlying rock volume (Figure 2). Despite the seafloor effect is obviously removed by the 4D approach, it was included in the single 3D models to evaluate the gravity signature related to bathymetry. Seafloor is found between 85 and 100 m below sea level, and the bathymetry shows a main step (∼10 m) in the eastern area of the modeled region and a gentle westward descending gradient across the entire area ( Figure 1B).
Within the geometric model, we set density values according to the petrophysical and geochemical data recovered from the reports available in the Volve database. In particular, by averaging the density values provided for the units surrounding Hugin, we set a bulk reference density for the volume of 2670 kg m ?3 . This value was used to calculate density contrasts of sea water volume (−1640 kg m −3 ) and the Hugin volume (−70 kg m −3 ) -i.e., we assume that seawater and Hugin have a density of 1030 and 2600 kg m −3 , respectively. The density contrasts were used as input for a pre-production 3D forward model. The procedure to calculate the forward gravity models of the volumes is based on the algorithm proposed by Li and Oldenburg (1998).
The vertical component of the gravity field produced by the density ρ(x,y,z) is given by: where V is the anomalous mass volume, r 0 = (x 0 , y 0 , z 0 ) is a vector locating the observation point, r = (x, y, z) locates the source, and γ is the gravitational constant.
Assuming a constant density contrast within each prismatic cell in a 3D orthogonal mesh, the gravity field at the ith observation location is given by: where ρ j and V j are the density and volume of the jth cell, respectively.
To evaluate density changes induced by fluid production and injection, we recovered data from the activity reports provided in the Volve dataset. The total reported volume balance accounting for 1.5 × 10 9 Sm 3 of volumes produced and 30 × 10 6 Sm 3 of injected water is 1.47 × 10 9 Sm 3 (Sm 3 stands for standard cubic meters and is referred to 15 • C and 1010 × 10 2 Pa) including oil, gas, and water production and injection. At reservoir conditions (106 • C and 3.28 × 10 7 Pa), this volume becomes 3.21 × 10 8 m 3 . At the same temperature and pressure conditions, the fluids in the reservoir have a density of 710 kg m −3 as reported in the geochemical analyses of samples from the discovery well 15/9-19-SR and included in the Volve dataset. Thus, the mass loss in the reservoir, accounting for production and injection activities and reservoir conditions, is of 2.28 × 10 11 kg. Considering that injection of water was contemporaneous to production (wells 15/9-F-4 and 15/9-F-5 were active between April 2008 and September 2016), we assume zero compaction of the Hugin sandstones during this period.
If a homogeneous mass loss across the entire reservoir is assumed, after the production each cell in the discretized Hugin volume should weigh 8.27 × 10 5 kg less than before productioni.e., the total mass loss divided by the total number of cells. In other words, knowing the volume of each cell (1.56 × 10 4 m 3 ), the production and injection activities resulted in a reservoir formation density loss of ∼53 kg m −3 . We note here that this configuration of the post-production model provides no preferential "paths" for density variations within the reservoir because it is not including constraints regarding injection and production wells' locations. Finally, the achieved density loss value is used to correct the pre-production density of the Hugin volume and perform a post-production 3D gravity forward model of the total volume that is compared with the preproduction model to evaluate gravity changes related to the production activity. Figure 3 shows the gravity signatures obtained after 3D forward calculation of the pre-production model (Figure 3a) and of the post-production model (Figure 3b) and the differential signature obtained by removing the forward-calculated gravity of the preproduction gravity from the post-production (Figure 3c). All these maps are produced with regular station spacing of 50 m located on sea surface, but similar results are obtained using a double station spacing of 100 m.

RESULTS
Despite the results of the pre-and post-production, forward models are extremely similar in the way that it would not be possible to qualitatively locate any local undulation; the difference between these two gravity data highlights a ∼−13 µGal (−130 × 10 −9 m s −2 ) minimum centered at 435256 E, 6478496 N (ED50, UTM 31N). Considering the modeled volume and the imposed density change, we interpret this minimum as representing the gravity effect of oil, gas, and water production and injection activities at the Volve field between February 2008 and September 2016.
The resulting residuals are in the same order of magnitude than those retrieved in the Sleipner field through gravity data acquisition (e.g., Alnes et al., 2008Alnes et al., , 2011. In a final processing step, we calculate the tilt derivative (Miller and Singh, 1994a;Verduzco et al., 2004) of the anomaly in Figure 3c. The horizontal gradient of the tilt derivative (THDR) is computed in order to locate zero values of the THDR representing the boundary of the causative source (e.g., Miller and Singh, 1994b). Projecting the calculated THDR (Figure 3d) on the reservoir top (Figure 4), we note that the THDR values ≈0 are surrounded by the most producing wells -i.e., 15/9-11, 15/9-12, and 15/9-14 whose production accounts for more than 98% of the total volume ( Figure 4B). This evidence supports the former interpretation concerning the gravity minima observed in Figure 3c being caused by reservoir production and injection activities.

DISCUSSION
The 4D gravity approach is not a new idea to monitor fluid-related density variations, but its application has always relied on high-precision gravity data acquisitions repeated in time and on source depth <2500 m (e.g., Ferguson et al., 2007;Vasilevskiy and Dashevsky, 2007;Alnes et al., 2008Alnes et al., , 2011Elliott and Braun, 2016). In this work, we applied a simple modeling procedure to evaluate the 4D gravity response computed from 3D forward models repeated at pre-production (February 2008) and post-production (September 2016) time steps. This procedure allowed to evaluate the gravity effects of pre-determined density changes in the >2750 m deep Volve reservoir without gravity data acquisition. The procedure was constrained by production data and geometry of the reservoir resulting from seismic and well data.
We assumed sealing conditions of the bounding faults and isotropy of the reservoir volume. Without detailed data about fault behavior and anisotropies within the reservoir levels, we made these assumptions to confine the mass variation within the modeled reservoir and to homogeneously distribute the mass change within the entire reservoir volume. However, the reference volume can be further parametrized in order to locate density changes according to available data or simulated scenarios. For example, to investigate eventual effects on the spatial distribution of density changes after introducing a leaking fault in the model, the reservoir volume should be extended beyond the investigated leaking fault and be parametrized with proper pre-production and post-production density contrasts. Furthermore, volume parametrization allows to model also eventual noise sources -e.g., topography/bathymetry, local density changes in surrounding volumes, as long as these are quantifiable by density or volume changes within the modeled period.
Based on the well-known algorithm proposed by Li and Oldenburg (1998), we modeled the gravity signature without filtering any regional or topographic components that are naturally removed by the differential approach. The only requirements are given by well-constrained geometries of the reservoir, availability of data about produced and injected volumes, and petrophysical properties of the fluids and reservoir formation. All these data are normally produced for oil and gas field characterization.
Interestingly, the procedure locates the causative source of the residual anomaly (Figure 3c), within the area of maximum production, where THDR values ≈ 0 are surrounded by the three most productive wells (Figure 4). In general, it is not surprising that a gravity low is centered over its causative source's top as we observed in the Volve field (Figure 4). However, what is surprising in this case is the observation per se of a minima being caused by ∼−53 kg m −3 density change at reservoir conditions and whose source is below 2750 m depth. In fact, the Volve field represents a challenging study case for the gravity method because, with a size/depth ratio slightly below the critical value of 1, gravity signatures related to density changes within the modeled volume should not be detectable because of the rapid attenuation of the gravity signal with increasing distance between the observation point and the source -e.g., Stenvold et al. (2008) and references therein.
This outcome strengthens the conclusions obtained by Stenvold et al. (2008) by expanding the retrievable depth for fluid-related density changes below the threshold of 2500 m. Furthermore, the Volve case modeling demonstrates that even without a detailed parametrization of the reservoir units, the simple forward calculation based on reservoir geometry and produced and injected volumes may prove reliable. The applicability of this method depends on the size/depth ratio of the source and on the magnitude of the density change both in terms of the density contrast imposed in the single cell and of thickness of the modeled reservoir carrying a density contrast. Approximating the reservoir shape to a horizontal slice of a cylindrical source (Telford et al., 1990;Stenvold et al., 2008) with radius 1.5 km, the imposed density change of 53 kg m −3 at a depth of 2750 m should produce a 10 µGal gravity change with reservoir thickness of ∼40 m (eq. 1 in Stenvold et al., 2008). Following the same approximation, a reservoir thickness of 100 m with the same density change and the same radius should produce a gravity change of ∼10 µGal at a depth of ∼4000 m. Of course, these approximations are not valid when dealing with complex structures and morphologies such as those often found in real geological contexts like the Volve field.
Considering the position of the gravity minimum and that the post-production input model assumed a homogeneous density change within the entire reservoir volume without constraints regarding injection and production wells, this outcome strengthens the reliability of the retrieved gravity anomaly. Moreover, we expect this approach to retrieve higher amplitude differential gravity residuals when dealing with higher density contrasts than those modeled at the Volve field. In fact, a scenario where reservoir fluids (typically oil, water, and gas) are replaced by CO 2 sequestration or gas storage should produce higher density contrasts and in turn should allow to resolve density changes related to gas propagation better.

CONCLUSION
In this work, we successfully locate the causative source for −53 kg m −3 density changes due to fluid production from a reservoir by 3D forward calculation of gravity at pre-and post-production time steps. The Volve field was chosen as test site because it represents a challenging reservoir with a size/depth ratio slightly below 1, and in these conditions, gravity methods should not allow detection of fluid-induced small density changes. Moreover, the Volve field represents a unique opportunity because in this case, data that usually are kept confidential also after exploitation have been made available to the community allowing the parametrization of the models.
The density change used as input for the post-production model, once removed the pre-production gravity signature, resulted in differential gravity residuals with a −13 µGal minimum above the Volve field. Unsurprisingly, despite the input density variation was distributed within the entire reservoir volume, THDR minimum located on top of the reservoir, between the three most productive wells in the field, whose total production accounted for more than 98% of total mass balance.
This work demonstrates that the 4D gravity method, whether it is constrained to gravity data acquisitions or to forward modeling of known geometries and density changes, represents a reliable technique even in challenging cases. Interestingly, the production-related anomalies observed at 2750 m depth in the Volve field should be theoretically retrievable by seafloor gravimeters with 5 µGal accuracy. However, this goal is still difficult to achieve and likely requires significant efforts in terms of number of stations in order to cover small sources like the Volve field.
Finally, relying on the availability of data that are confidential in the case of oil fields or very expensive to produce for research purposes, this approach has intrinsic limitations in the possibility of being tested in more complex geological scenarios. However, the owners of such data should follow the example of Equinor and release data regarding fields that are no longer exploited; this will likely open new ways in understanding geological and geophysical processes.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. These data can be found here: Volve data village page: https://www. equinor.com/en/how-and-why/digitalisation-in-our-dna/volvefield-data-village-download.html.

AUTHOR CONTRIBUTIONS
PM conceived the study, produced the models, and wrote the manuscript.

FUNDING
This work was supported by the Chieti-Pescara University funds to PM.