Well Log Fluid Substitution¶

In [1]:

import numpy as np
import pandas as pd
import os
import pathlib
from digirock.datasets import fetch_example_data

import matplotlib.pyplot as plt

import numpy as np import pandas as pd import os import pathlib from digirock.datasets import fetch_example_data import matplotlib.pyplot as plt

Fetch and setup the example data, which includes well F-4 from the Volve field.

In [2]:

data = fetch_example_data()

data = fetch_example_data()

In [3]:

logs = pd.read_csv(data["volve_f4.csv"], index_col=0)
logs["VELP"] = 304800.0 / logs["DT"] # get to meters per second
logs["VELS"] = 304800.0 / logs["DTS"] # get to meters per second

# load some formation volume factor tables for the oil
bo_tab = pd.read_csv(data["PVT_BO.inc"], delim_whitespace=True, header=None, names=["pres", "bo"])
rs_tab = pd.read_csv(data["PVT_RS.inc"], delim_whitespace=True, header=None, names=["pres", "rs"])

# pres is in bar -> convert to MPa for digirock
bo_tab["pres"] = bo_tab["pres"] * 0.1
rs_tab["pres"] = rs_tab["pres"] * 0.1

logs = pd.read_csv(data["volve_f4.csv"], index_col=0) logs["VELP"] = 304800.0 / logs["DT"] # get to meters per second logs["VELS"] = 304800.0 / logs["DTS"] # get to meters per second # load some formation volume factor tables for the oil bo_tab = pd.read_csv(data["PVT_BO.inc"], delim_whitespace=True, header=None, names=["pres", "bo"]) rs_tab = pd.read_csv(data["PVT_RS.inc"], delim_whitespace=True, header=None, names=["pres", "rs"]) # pres is in bar -> convert to MPa for digirock bo_tab["pres"] = bo_tab["pres"] * 0.1 rs_tab["pres"] = rs_tab["pres"] * 0.1

In [4]:

logs.describe()

logs.describe()

Out[4]:

	DT	DTS	TVDSS	GR	NPHI	DENS	XLOC	YLOC	SW_SURVEY1	SW_SURVEY2	PERMZ	PORO	VSH	OWT	VELP	VELS
count	4815.000000	4246.000000	5988.000000	5988.000000	4926.000000	4964.000000	5988.000000	5.988000e+03	1138.000000	1138.000000	1070.000000	1070.000000	5988.000000	5988.000000	4815.000000	4246.000000
mean	78.140018	142.814872	2771.805207	46.109442	0.180302	2.441509	435690.865810	6.477987e+06	0.295872	0.646379	548.386215	0.229878	0.260368	1262.285207	4008.790939	2159.695890
std	13.415549	15.671214	200.146933	44.143014	0.097218	0.166057	33.141649	1.578144e+02	0.295049	0.259817	531.485748	0.025175	0.217152	57.942543	644.440103	233.760768
min	53.358002	108.222000	2400.089355	1.000000	0.045700	1.334100	435601.125000	6.477678e+06	0.042300	0.142800	0.491576	0.130000	0.020500	1140.923828	2331.025434	1601.995596
25%	67.192200	130.147629	2599.324951	15.153900	0.107500	2.300000	435670.031250	6.477854e+06	0.131600	0.684475	62.688911	0.216640	0.080800	1217.474945	3499.622661	1978.556591
50%	76.447197	142.112801	2788.841186	35.016998	0.157000	2.496600	435699.500000	6.478012e+06	0.167800	0.764000	513.006897	0.230275	0.209550	1262.712891	3987.065741	2144.775166
75%	87.095104	154.051697	2951.317078	56.746299	0.227675	2.560025	435719.468750	6.478132e+06	0.248900	0.779200	813.894226	0.248776	0.359550	1315.449036	4536.240833	2341.955849
max	130.757904	190.262695	3083.119385	469.738007	0.732400	3.201300	435727.593750	6.478198e+06	1.000000	1.000000	2000.000000	0.281082	0.799800	1350.738525	5712.357851	2816.432885

Setup the PEM¶

We need to first build the basic blocks of our rock model including an oil/water WoodsFluid and Voigt-Reuss-Hill grain matrix.

In [5]:

from digirock import WaterBW92, OilPVT, WoodsFluid, Mineral, VRHAvg

# FLUID PROPS
temp = 110  # reservoir temperature in degC
sal = 151200  # reservoir brine salinity ppm
pres = 32 # reseroir insitu pressure MPa

volve_brine = WaterBW92(salinity=sal)
volve_oil = OilPVT(std_density=0.885, gas_sg=1.09956)
volve_oil.set_pvt(140, bo_tab["bo"].values, bo_tab["pres"].values)

fl = WoodsFluid(["sw", "so"], [volve_brine, volve_oil])

# MATRIX PROPS
volve_sand = Mineral(2.634, 31.55, 29.69, 'sand')
volve_clay = Mineral(2.78, 28.45, 7.95, 'clay')

vrh_avg = VRHAvg(["vsand", "vclay"], [volve_sand, volve_clay])

from digirock import WaterBW92, OilPVT, WoodsFluid, Mineral, VRHAvg # FLUID PROPS temp = 110 # reservoir temperature in degC sal = 151200 # reservoir brine salinity ppm pres = 32 # reseroir insitu pressure MPa volve_brine = WaterBW92(salinity=sal) volve_oil = OilPVT(std_density=0.885, gas_sg=1.09956) volve_oil.set_pvt(140, bo_tab["bo"].values, bo_tab["pres"].values) fl = WoodsFluid(["sw", "so"], [volve_brine, volve_oil]) # MATRIX PROPS volve_sand = Mineral(2.634, 31.55, 29.69, 'sand') volve_clay = Mineral(2.78, 28.45, 7.95, 'clay') vrh_avg = VRHAvg(["vsand", "vclay"], [volve_sand, volve_clay])

Create the Gassmann substitution Blender¶

To perform the fluid substitution, a class which blends the props data is needed. Here we build a simple example of a GassmannSub class but it could be arbitrarilly more complex.

In [6]:

from digirock import Blend
from digirock.models import dryframe_acoustic, gassmann_fluidsub
from digirock.elastic import acoustic_bulk_moduli, acoustic_shear_moduli, acoustic_velp, acoustic_vels

class GassmannSub(Blend):
    
    _methods = ["vp", "vs", "density", "bulk_modulus", "shear_modulus"]
    
    def __init__(
        self,
        zero_porosity_model,
        fluid_model,
        sw_key="sw",
        vp_key="vp",
        vs_key="vs",
        rho_key="rho",
        poro_key="poro", 
        name=None
    ):
        elements = [zero_porosity_model, fluid_model]
        blend_keys = [sw_key, vp_key, vs_key, rho_key, poro_key]
        super().__init__(blend_keys, elements, self._methods, name=name)
        
    def _get_elastic_keys(self):
        return self.blend_keys[1:4]
    
    def dry_bulk_modulus(self, props, **kwargs):
        kfl = self.elements[1].bulk_modulus(props, **kwargs)
        k0 = self.elements[0].bulk_modulus(props, **kwargs)
        vpk, vsk, rk = self._get_elastic_keys()
        porok = self.blend_keys[4]
        ksat = acoustic_bulk_moduli(props[vpk], props[vsk], props[rk])
        kdry = dryframe_acoustic(ksat, kfl, k0, props[porok])
        return kdry
    
    def bulk_modulus(self, props, props2, **kwargs):
        kdry = self.dry_bulk_modulus(props)
        kfl = self.elements[1].bulk_modulus(props2, **kwargs)
        k0 = self.elements[0].bulk_modulus(props2, **kwargs)
        porok = self.blend_keys[4]
        return gassmann_fluidsub(kdry, kfl, k0, props2[porok])
    
    def shear_modulus(self, props, **kwargs):
        return acoustic_shear_moduli(props[self.blend_keys[2]], props[self.blend_keys[3]])
    
    def density(self, props, props2, **kwargs):
        porok = self.blend_keys[4]
        vpk, vsk, rk = self._get_elastic_keys()
        kfl1 = self.elements[1].density(props, **kwargs)
        den_min = props[rk] - props[porok] * kfl1
        kfl2 = self.elements[1].density(props2, **kwargs)
        return den_min + props2[porok] * kfl2
    
    def vp(self, props, props2, **kwargs):
        k = self.bulk_modulus(props, props2, **kwargs)
        mu = self.shear_modulus(props, **kwargs)
        den = self.density(props, props2, **kwargs)
        return acoustic_velp(k, mu, den)

    def vs(self, props, props2, **kwargs):
        mu = self.shear_modulus(props, **kwargs)
        den = self.density(props, props2, **kwargs)
        return acoustic_vels(mu, den)
        
    
test_props = {"sw":0.5, "pres":pres, "temp":temp, "vclay":0.5, "poro":0.2, "vp":3500, "vs":2500, "rho":2.2}
test_props2 = {"sw":1.0, "pres":pres, "temp":temp, "vclay":0.5, "poro":0.2}

gassub = GassmannSub(vrh_avg, fl)
print("Dry Bulk Mod: ", gassub.dry_bulk_modulus(test_props))
print("Subbed Bulk Mod: ", gassub.bulk_modulus(test_props, test_props2))
print("Subbed Shear Mod: ", gassub.shear_modulus(test_props))
print("Subbed Density: ", gassub.density(test_props, test_props2))
print("Subbed VP: ", gassub.vp(test_props, test_props2))
print("Subbed VS: ", gassub.vs(test_props, test_props2))

from digirock import Blend from digirock.models import dryframe_acoustic, gassmann_fluidsub from digirock.elastic import acoustic_bulk_moduli, acoustic_shear_moduli, acoustic_velp, acoustic_vels class GassmannSub(Blend): _methods = ["vp", "vs", "density", "bulk_modulus", "shear_modulus"] def __init__( self, zero_porosity_model, fluid_model, sw_key="sw", vp_key="vp", vs_key="vs", rho_key="rho", poro_key="poro", name=None ): elements = [zero_porosity_model, fluid_model] blend_keys = [sw_key, vp_key, vs_key, rho_key, poro_key] super().__init__(blend_keys, elements, self._methods, name=name) def _get_elastic_keys(self): return self.blend_keys[1:4] def dry_bulk_modulus(self, props, **kwargs): kfl = self.elements[1].bulk_modulus(props, **kwargs) k0 = self.elements[0].bulk_modulus(props, **kwargs) vpk, vsk, rk = self._get_elastic_keys() porok = self.blend_keys[4] ksat = acoustic_bulk_moduli(props[vpk], props[vsk], props[rk]) kdry = dryframe_acoustic(ksat, kfl, k0, props[porok]) return kdry def bulk_modulus(self, props, props2, **kwargs): kdry = self.dry_bulk_modulus(props) kfl = self.elements[1].bulk_modulus(props2, **kwargs) k0 = self.elements[0].bulk_modulus(props2, **kwargs) porok = self.blend_keys[4] return gassmann_fluidsub(kdry, kfl, k0, props2[porok]) def shear_modulus(self, props, **kwargs): return acoustic_shear_moduli(props[self.blend_keys[2]], props[self.blend_keys[3]]) def density(self, props, props2, **kwargs): porok = self.blend_keys[4] vpk, vsk, rk = self._get_elastic_keys() kfl1 = self.elements[1].density(props, **kwargs) den_min = props[rk] - props[porok] * kfl1 kfl2 = self.elements[1].density(props2, **kwargs) return den_min + props2[porok] * kfl2 def vp(self, props, props2, **kwargs): k = self.bulk_modulus(props, props2, **kwargs) mu = self.shear_modulus(props, **kwargs) den = self.density(props, props2, **kwargs) return acoustic_velp(k, mu, den) def vs(self, props, props2, **kwargs): mu = self.shear_modulus(props, **kwargs) den = self.density(props, props2, **kwargs) return acoustic_vels(mu, den) test_props = {"sw":0.5, "pres":pres, "temp":temp, "vclay":0.5, "poro":0.2, "vp":3500, "vs":2500, "rho":2.2} test_props2 = {"sw":1.0, "pres":pres, "temp":temp, "vclay":0.5, "poro":0.2} gassub = GassmannSub(vrh_avg, fl) print("Dry Bulk Mod: ", gassub.dry_bulk_modulus(test_props)) print("Subbed Bulk Mod: ", gassub.bulk_modulus(test_props, test_props2)) print("Subbed Shear Mod: ", gassub.shear_modulus(test_props)) print("Subbed Density: ", gassub.density(test_props, test_props2)) print("Subbed VP: ", gassub.vp(test_props, test_props2)) print("Subbed VS: ", gassub.vs(test_props, test_props2))

Dry Bulk Mod:  [6.15444211]
Subbed Bulk Mod:  [9.21680372]
Subbed Shear Mod:  13.750000000000002
Subbed Density:  [2.23761226]
Subbed VP:  [3508.88766769]
Subbed VS:  [2478.89957096]

In [7]:

gassub.tree

gassub.tree

Out[7]:

GassmannSub_140699708089200
├── class : <class '__main__.GassmannSub'>
├── props_keys : ['sw', 'vp', 'vs', 'rho', 'poro']
├── blend_keys : ['sw', 'vp', 'vs', 'rho', 'poro']
├── methods : ('vp', 'vs', 'density', 'bulk_modulus', 'shear_modulus')
├── n_elements : 2
├── VRHAvg_140699708275440
│   ├── class : <class 'digirock._frames._frames.VRHAvg'>
│   ├── props_keys : ['vsand', 'vclay']
│   ├── blend_keys : ['vsand', 'vclay']
│   ├── methods : ('density', 'vp', 'vs', 'shear_modulus', 'bulk_modulus')
│   ├── n_elements : 2
│   └── elements
│       ├── sand
│       │   ├── class : <class 'digirock._frames._minerals.Mineral'>
│       │   ├── props_keys : []
│       │   ├── bulk_modulus : 31.55
│       │   ├── shear_modulus : 29.69
│       │   ├── dens : 2.634
│       │   ├── vp : 5196.834306908865
│       │   └── vs : 3357.3546009435527
│       └── clay
│           ├── class : <class 'digirock._frames._minerals.Mineral'>
│           ├── props_keys : []
│           ├── bulk_modulus : 28.45
│           ├── shear_modulus : 7.95
│           ├── dens : 2.78
│           ├── vp : 3747.9010912680255
│           └── vs : 1691.0683694681975
└── WoodsFluid_140699708274480
    ├── class : <class 'digirock._fluids._fluid_mod.WoodsFluid'>
    ├── props_keys : ['sw', 'so']
    ├── blend_keys : ['sw', 'so']
    ├── methods : ('density', 'bulk_modulus', 'shear_modulus', 'velocity')
    ├── n_elements : 2
    └── elements
        ├── WaterBW92_140699708272752
        │   ├── class : <class 'digirock._fluids._water.WaterBW92'>
        │   ├── props_keys : ['temp', 'pres']
        │   └── salinity : 0.1512
        └── OilPVT_140699708274528
            ├── class : <class 'digirock._fluids._oil.OilPVT'>
            ├── props_keys : ['bo']
            ├── api : 28.38700564971751
            ├── std_density : 0.885
            ├── gas_sg : 1.09956
            └── pvt
                ├── type : fixed_rs
                ├── rs : 140
                ├── bo : table
                ├── pres : table
                ├── bo_table : <xarray.DataArray (rs: 1, pres: 14)>
                │   array([[1.17 , 1.248, 1.322, 1.398, 1.483, 1.477, 1.467, 1.456, 1.445,
                │           1.435, 1.427, 1.419, 1.411, 1.405]])
                │   Coordinates:
                │     * rs       (rs) int64 140
                │     * pres     (pres) float64 5.0 10.0 15.0 20.0 25.0 ... 45.0 50.0 55.0 
                │   60.0 65.0
                ├── bo_min : 1.17
                ├── bo_max : 1.483
                ├── pres_min : 5.0
                └── pres_max : 65.0

Performing fluid substitution on the logs¶

We simply need to define the insitu state props1 and the state at which we wish to calculate the new logs props2. And these can be passed to gassub, our model class, which handles the details of the fluid substitution.

In [8]:

fig, axs = plt.subplots(ncols=5, figsize=(10,6), sharey=True)

sub = logs[np.logical_and(logs.index > 3250, logs.index < 3420)]

axs[0].plot(sub["VSH"], -sub.index, color="green")
axs[1].plot(sub["SW_SURVEY1"], -sub.index)
axs[1].plot(sub["SW_SURVEY2"], -sub.index)
axs[2].plot(sub["VELP"], -sub.index)
axs[3].plot(sub["VELS"], -sub.index)
axs[4].plot(sub["DENS"], -sub.index)

props1 = {
    "pres":pres, 
    "temp":temp, 
    "poro":sub["PORO"].values, 
    "vclay":sub["VSH"].values,
    "vp":sub["VELP"].values,
    "vs":sub["VELS"].values,
    "rho":sub["DENS"].values,
    "sw":sub["SW_SURVEY1"].values
}
props2 = props1.copy()
props2["sw"] = sub["SW_SURVEY2"].values

vp2 = gassub.vp(props1, props2)
vs2 = gassub.vs(props1, props2)
rho2 = gassub.density(props1, props2)

axs[2].plot(vp2, -sub.index)
axs[3].plot(vs2, -sub.index)
axs[4].plot(rho2, -sub.index)

axs[0].set_title("VSH")
axs[1].set_title("SW1, SW2")
axs[2].set_title("VP1, VP2")
axs[3].set_title("VS1, VS2")
axs[4].set_title("RHOB1, RHOB2")

fig, axs = plt.subplots(ncols=5, figsize=(10,6), sharey=True) sub = logs[np.logical_and(logs.index > 3250, logs.index < 3420)] axs[0].plot(sub["VSH"], -sub.index, color="green") axs[1].plot(sub["SW_SURVEY1"], -sub.index) axs[1].plot(sub["SW_SURVEY2"], -sub.index) axs[2].plot(sub["VELP"], -sub.index) axs[3].plot(sub["VELS"], -sub.index) axs[4].plot(sub["DENS"], -sub.index) props1 = { "pres":pres, "temp":temp, "poro":sub["PORO"].values, "vclay":sub["VSH"].values, "vp":sub["VELP"].values, "vs":sub["VELS"].values, "rho":sub["DENS"].values, "sw":sub["SW_SURVEY1"].values } props2 = props1.copy() props2["sw"] = sub["SW_SURVEY2"].values vp2 = gassub.vp(props1, props2) vs2 = gassub.vs(props1, props2) rho2 = gassub.density(props1, props2) axs[2].plot(vp2, -sub.index) axs[3].plot(vs2, -sub.index) axs[4].plot(rho2, -sub.index) axs[0].set_title("VSH") axs[1].set_title("SW1, SW2") axs[2].set_title("VP1, VP2") axs[3].set_title("VS1, VS2") axs[4].set_title("RHOB1, RHOB2")

Out[8]:

Text(0.5, 1.0, 'RHOB1, RHOB2')