# framework
using GeoStats
# IO modules
using GeoIO
# viz modules
import CairoMakie as Mke14 Petroleum reservoirs
Petroleum reservoirs present various modeling challenges related to their complex geometry and distribution of rock and fluid properties. Some of these challenges are still open in industry due the lack of software for advanced geospatial data science with unstructured meshes. In this chapter, we will illustrate how an important “oil-in-place” calculation in reservoir management can be automated with the framework.
TOOLS COVERED: @groupby, @transform, @combine, Unitify, Unit, GHC, volume, viewer
MODULES:
Note
Although we use CairoMakie.jl in this book, many of the 3D visualizations in this chapter demand a more performant Makie.jl backend. Consider using GLMakie.jl if you plan to reproduce the code locally.
14.1 Data
We will use reservoir simulation results of the Norne benchmark case, a real oil field from the Norwegian Sea. For more information, please check the OPM project. These results were simulated with the JutulDarcy.jl reservoir simulator by Møyner (2024).
In particular, we will consider only two time steps of the simulation, named norne1.vtu and norne2.vtu. The data are stored in the open VTK format with .vtu extension, indicating that it is georeferenced over an unstructured mesh:
norne₁ = GeoIO.load("data/norne1.vtu")
norne₂ = GeoIO.load("data/norne2.vtu")
norne₁ |> viewer
Note
The vtk.org website provides official documentation for the various VTK file formats, including formats for image data (.vti), rectilinear grids (.vtr), structured grids (.vts), unstructured meshes (.vtu), etc.
14.2 Objectives
The Volume of Oil In Place (\(VOIP\)) is a global estimate of the volume of oil trapped in the subsurface. It is defined as an integral over the volume \(V\) of the reservoir:
\[ VOIP = \int_V S_o \phi\ dV \]
where \(\phi\) is the rock porosity and \(S_o\) is the oil saturation. The integrand can be converted into Mass of Oil in Place (\(MOIP\)) using the oil density \(\rho_o\):
\[ MOIP = \int_V \rho_o S_o \phi\ dV \]
Likewise, the Mass of Water In Place (\(MWIP\)) and Mass of Gas in Place (\(MGIP\)) are defined using the respective fluid saturations and densities.
Our main objective is to estimate these masses of fluids in place over a reservoir model with non-trivial geometry, for different time steps within a physical reservoir simulation. This can be useful to understand rates of depletion and guide reservoir management.
Secondary objectives include the localization (through 3D visualization) of zones with high mass of hydrocarbons (oil + gas), and the calculation of zonal depletion, i.e., the change of hydrocarbon mass per zone, from a time step \(t_1\) to a time step \(t_2\):
\[ Depletion = \left\{MOIP + MGIP\right\}_{t_1} - \left\{MOIP + MGIP\right\}_{t_2} \]
14.3 Methodology
In order to identify zones of the reservoir with high mass of hydrocarbons, we need to compute the fluids in place for each element of the mesh, and group the elements based on their calculated masses. Given the zones, we can compute the zonal depletion.
The proposed methodology has the following steps:
- Analysis of oil, gas and water in place
- Localization of hydrocarbon zones
- Calculation of zonal depletion
14.3.1 Fluid analysis
Before we start our calculations, we need to rename the variables in the dataset to match our concise notation. We also need to correct the units of the variables to make sure that our final report has values that are easy to read.
The following pipeline performs the desired cleaning steps by exploiting bracket notation (e.g., [kg/m^3]) for units. The Unitify transform takes a geotable with bracket notation as input and converts the values of columns to unitful values:
clean = Select(
"porosity" => "ϕ",
"saturation_oil" => "So",
"saturation_gas" => "Sg",
"saturation_water" => "Sw",
"density_oil" => "ρo [kg/m^3]",
"density_gas" => "ρg [kg/m^3]",
"density_water" => "ρw [kg/m^3]"
) → Unitify()SequentialTransform
├─ Select(selector: [:porosity, :saturation_oil, :saturation_gas, :saturation_water, :density_oil, :density_gas, :density_water], newnames: [:ϕ, :So, :Sg, :Sw, Symbol("ρo [kg/m^3]"), Symbol("ρg [kg/m^3]"), Symbol("ρw [kg/m^3]")])
└─ Unitify()The resulting geotable has variables with concise names and correct units:
reservoir₁ = clean(norne₁)
reservoir₂ = clean(norne₂)| ϕ | So | Sg | Sw | ρo | ρg | ρw | geometry |
|---|---|---|---|---|---|---|---|
| Continuous | Continuous | Continuous | Continuous | Continuous | Continuous | Continuous | Hexahedron |
| [NoUnits] | [NoUnits] | [NoUnits] | [NoUnits] | [kg m^-3] | [kg m^-3] | [kg m^-3] | 🖈 Cartesian{NoDatum} |
| 0.30081 | 0.894952 | 0.0 | 0.105049 | 718.387 kg m^-3 | 189.051 kg m^-3 | 994.545 kg m^-3 | Hexahedron((x: 455586.0 m, y: 7.32115e6 m, z: 2586.74 m), ..., (x: 455636.0 m, y: 7.32122e6 m, z: 2605.72 m)) |
| 0.299529 | 0.000697226 | 0.894094 | 0.105209 | 716.681 kg m^-3 | 188.769 kg m^-3 | 994.486 kg m^-3 | Hexahedron((x: 455608.0 m, y: 7.32107e6 m, z: 2568.85 m), ..., (x: 455659.0 m, y: 7.32114e6 m, z: 2585.64 m)) |
| 0.29655 | 0.00040746 | 0.894254 | 0.105339 | 716.723 kg m^-3 | 188.685 kg m^-3 | 994.457 kg m^-3 | Hexahedron((x: 455625.0 m, y: 7.32099e6 m, z: 2563.06 m), ..., (x: 455678.0 m, y: 7.32106e6 m, z: 2573.92 m)) |
| 0.29324 | 0.000401796 | 0.894189 | 0.10541 | 716.745 kg m^-3 | 188.64 kg m^-3 | 994.442 kg m^-3 | Hexahedron((x: 455650.0 m, y: 7.32092e6 m, z: 2559.15 m), ..., (x: 455706.0 m, y: 7.32099e6 m, z: 2569.36 m)) |
| 0.289647 | 0.00040205 | 0.894385 | 0.105212 | 716.765 kg m^-3 | 188.6 kg m^-3 | 994.427 kg m^-3 | Hexahedron((x: 455677.0 m, y: 7.32085e6 m, z: 2555.98 m), ..., (x: 455731.0 m, y: 7.32091e6 m, z: 2567.17 m)) |
| 0.28579 | 0.000402928 | 0.894402 | 0.105195 | 716.783 kg m^-3 | 188.565 kg m^-3 | 994.415 kg m^-3 | Hexahedron((x: 455700.0 m, y: 7.32078e6 m, z: 2553.79 m), ..., (x: 455756.0 m, y: 7.32084e6 m, z: 2562.11 m)) |
| 0.280003 | 0.000401767 | 0.894401 | 0.105197 | 716.788 kg m^-3 | 188.554 kg m^-3 | 994.411 kg m^-3 | Hexahedron((x: 455724.0 m, y: 7.3207e6 m, z: 2552.27 m), ..., (x: 455781.0 m, y: 7.32076e6 m, z: 2560.06 m)) |
| 0.274279 | 0.000383991 | 0.859452 | 0.140164 | 716.793 kg m^-3 | 188.546 kg m^-3 | 994.408 kg m^-3 | Hexahedron((x: 455751.0 m, y: 7.32063e6 m, z: 2551.07 m), ..., (x: 455804.0 m, y: 7.32068e6 m, z: 2561.31 m)) |
| 0.270403 | 0.000385463 | 0.859456 | 0.140159 | 716.805 kg m^-3 | 188.521 kg m^-3 | 994.399 kg m^-3 | Hexahedron((x: 455785.0 m, y: 7.32057e6 m, z: 2549.52 m), ..., (x: 455831.0 m, y: 7.32061e6 m, z: 2559.75 m)) |
| 0.266089 | 0.000382203 | 0.854474 | 0.145144 | 716.818 kg m^-3 | 188.494 kg m^-3 | 994.389 kg m^-3 | Hexahedron((x: 455823.0 m, y: 7.3205e6 m, z: 2547.65 m), ..., (x: 455866.0 m, y: 7.32054e6 m, z: 2556.82 m)) |
| ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
We @transform the reservoir and compute masses of fluids for each element of the mesh using the formulae in the beginning of the chapter:
mass(reservoir) = @transform(reservoir,
:MOIP = :ρo * :So * :ϕ * volume(:geometry),
:MGIP = :ρg * :Sg * :ϕ * volume(:geometry),
:MWIP = :ρw * :Sw * :ϕ * volume(:geometry)
)
mass₁ = mass(reservoir₁)
mass₂ = mass(reservoir₂)
mass₁ |> Select("MWIP") |> viewer
14.3.2 Hydrocarbon zones
We compute the mass of hydrocarbon in place \(MHIP\) as the sum of oil and gas in the first time step, and cluster it with geostatistical hierarchical clustering (GHC) (Fouedjio 2016). The method requires an approximate number of clusters that we set to \(k=3\) (low, medium and high values) and a maximum range of geospatial association that we set to \(\lambda = 500m\).
Additionally, we set an upper bound nmax=1000 on the number of elements used in the dissimilarity matrix computation and the option as="zone" to name the column with clustering results.
zones = @transform(mass₁, :MHIP = :MOIP + :MGIP) |>
Select("MHIP") |> GHC(3, 500u"m", nmax=1000, as="zone")
zones |> viewer
14.3.3 Zonal depletion
We concatenate all variables of interest in a single geotable to be able to use the geospatial split-apply-combine pattern, and compute the final summary table with statistics per zone:
carbon₁ = mass₁ |> Select("MOIP" => "MOIP₁", "MGIP" => "MGIP₁")
carbon₂ = mass₂ |> Select("MOIP" => "MOIP₂", "MGIP" => "MGIP₂")
data = [carbon₁ carbon₂ zones]| MOIP₁ | MGIP₁ | MOIP₂ | MGIP₂ | zone | geometry |
|---|---|---|---|---|---|
| Continuous | Continuous | Continuous | Continuous | Categorical | Hexahedron |
| [kg] | [kg] | [kg] | [kg] | [NoUnits] | 🖈 Cartesian{NoDatum} |
| 7.57425e6 kg | 0.0 kg | 7.56895e6 kg | 0.0 kg | 3 | Hexahedron((x: 455586.0 m, y: 7.32115e6 m, z: 2586.74 m), ..., (x: 455636.0 m, y: 7.32122e6 m, z: 2605.72 m)) |
| 0.0 kg | 2.43976e6 kg | 7141.61 kg | 2.41219e6 kg | 3 | Hexahedron((x: 455608.0 m, y: 7.32107e6 m, z: 2568.85 m), ..., (x: 455659.0 m, y: 7.32114e6 m, z: 2585.64 m)) |
| 0.0 kg | 2.24126e6 kg | 3835.96 kg | 2.21634e6 kg | 3 | Hexahedron((x: 455625.0 m, y: 7.32099e6 m, z: 2563.06 m), ..., (x: 455678.0 m, y: 7.32106e6 m, z: 2573.92 m)) |
| 0.0 kg | 2.16366e6 kg | 3652.63 kg | 2.13943e6 kg | 3 | Hexahedron((x: 455650.0 m, y: 7.32092e6 m, z: 2559.15 m), ..., (x: 455706.0 m, y: 7.32099e6 m, z: 2569.36 m)) |
| 0.0 kg | 2.21737e6 kg | 3746.6 kg | 2.19304e6 kg | 3 | Hexahedron((x: 455677.0 m, y: 7.32085e6 m, z: 2555.98 m), ..., (x: 455731.0 m, y: 7.32091e6 m, z: 2567.17 m)) |
| 0.0 kg | 2.25638e6 kg | 3821.64 kg | 2.23166e6 kg | 3 | Hexahedron((x: 455700.0 m, y: 7.32078e6 m, z: 2553.79 m), ..., (x: 455756.0 m, y: 7.32084e6 m, z: 2562.11 m)) |
| 0.0 kg | 1.91029e6 kg | 3226.36 kg | 1.88936e6 kg | 3 | Hexahedron((x: 455724.0 m, y: 7.3207e6 m, z: 2552.27 m), ..., (x: 455781.0 m, y: 7.32076e6 m, z: 2560.06 m)) |
| 0.0 kg | 1.78157e6 kg | 2993.03 kg | 1.76212e6 kg | 3 | Hexahedron((x: 455751.0 m, y: 7.32063e6 m, z: 2551.07 m), ..., (x: 455804.0 m, y: 7.32068e6 m, z: 2561.31 m)) |
| 0.0 kg | 1.70066e6 kg | 2868.5 kg | 1.68211e6 kg | 3 | Hexahedron((x: 455785.0 m, y: 7.32057e6 m, z: 2549.52 m), ..., (x: 455831.0 m, y: 7.32061e6 m, z: 2559.75 m)) |
| 0.0 kg | 1.57036e6 kg | 2642.13 kg | 1.55327e6 kg | 3 | Hexahedron((x: 455823.0 m, y: 7.3205e6 m, z: 2547.65 m), ..., (x: 455866.0 m, y: 7.32054e6 m, z: 2556.82 m)) |
| ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
The depletion per zone can be computed with
summary = @chain data begin
@groupby(:zone)
@transform(:delta = :MOIP₁ + :MGIP₁ - :MOIP₂ - :MGIP₂)
@combine(:depletion = sum(:delta))
end| zone | depletion | geometry |
|---|---|---|
| Categorical | Continuous | MultiPolyhedron |
| [NoUnits] | [kg] | 🖈 Cartesian{NoDatum} |
| 3 | 4.23182e10 kg | Multi(21205×Hexahedron) |
| 2 | 1.86485e10 kg | Multi(6463×Hexahedron) |
| 1 | 2.04584e9 kg | Multi(16763×Hexahedron) |
or in \(Mg\) (ton) after a change of Unit:
summary |> Unit("depletion" => u"Mg")| zone | depletion | geometry |
|---|---|---|
| Categorical | Continuous | MultiPolyhedron |
| [NoUnits] | [Mg] | 🖈 Cartesian{NoDatum} |
| 3 | 4.23182e7 Mg | Multi(21205×Hexahedron) |
| 2 | 1.86485e7 Mg | Multi(6463×Hexahedron) |
| 1 | 2.04584e6 Mg | Multi(16763×Hexahedron) |
14.4 Summary
In this chapter, we illustrated the application of the framework in the petroleum industry. Among other things, we learned how to
- Perform simple calculations involving fluids in place and unstructured meshes.
- Identify zones of a petroleum reservoir using clustering methods and visualizations.
This open source technology can be used to create advanced dashboards for reservoir management without advanced programming skills. It addresses real issues raised by geospatial data scientists in industry who feel unproductive using rigid geomodeling software.