Frames¶

RockFrame classes are used in digirock to build the dry rock frames. They inherit the Blend class.

In [1]:

from digirock.utils import create_orthogonal_props
import numpy as np

%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib import rc

rc("font", size=14)
figsize = (15, 5)

# a function for repeat plotting at stages of development
def plot_elastic(element, props):
    fig, axes = plt.subplots(ncols=3, figsize=(8, 2.5), sharex=True)
    # Plot density
    axes[0].plot(props["VCLAY"], element.density(props))
    axes[0].set_ylabel("Density (g/cc)")
    # Plot moduli
    axes[1].plot(props["VCLAY"], element.bulk_modulus(props), label="bulk")
    axes[1].set_ylabel("Modulus (GPa)")
    axes[1].plot(props["VCLAY"], element.shear_modulus(props), label="shear")
    axes[1].legend()
    # Plot velocity
    axes[2].plot(props["VCLAY"], element.vp(props), label="vp")
    axes[2].set_ylabel("Velocity (m/s)")
    axes[2].plot(props["VCLAY"], element.vs(props), label="vs")
    axes[2].legend()

    for ax in axes:
        ax.set_xlabel("VCLAY (frac)")
    fig.tight_layout()

from digirock.utils import create_orthogonal_props import numpy as np %matplotlib inline import matplotlib.pyplot as plt from matplotlib import rc rc("font", size=14) figsize = (15, 5) # a function for repeat plotting at stages of development def plot_elastic(element, props): fig, axes = plt.subplots(ncols=3, figsize=(8, 2.5), sharex=True) # Plot density axes[0].plot(props["VCLAY"], element.density(props)) axes[0].set_ylabel("Density (g/cc)") # Plot moduli axes[1].plot(props["VCLAY"], element.bulk_modulus(props), label="bulk") axes[1].set_ylabel("Modulus (GPa)") axes[1].plot(props["VCLAY"], element.shear_modulus(props), label="shear") axes[1].legend() # Plot velocity axes[2].plot(props["VCLAY"], element.vp(props), label="vp") axes[2].set_ylabel("Velocity (m/s)") axes[2].plot(props["VCLAY"], element.vs(props), label="vs") axes[2].legend() for ax in axes: ax.set_xlabel("VCLAY (frac)") fig.tight_layout()

In [2]:

from digirock import VRHAvg, FixedPoroAdjust, NurCriticalPoroAdjust, LGStressModel, MacBethStressAdjust, Mineral

from digirock import VRHAvg, FixedPoroAdjust, NurCriticalPoroAdjust, LGStressModel, MacBethStressAdjust, Mineral

Create two end point minerals for a our sand clay system.

In [3]:

sand = Mineral(2.7, 29.9, 22.5)
clay = Mineral(2.51, 22, 9.4)

sand = Mineral(2.7, 29.9, 22.5) clay = Mineral(2.51, 22, 9.4)

Create a VRHAvg mixing model.

In [4]:

# the model -> rh
rf = VRHAvg(["VSAND", "VCLAY"], [sand, clay])

# check modulus variation with volume fraction
props = dict(VCLAY=np.linspace(0, 1, 51),)
plot_elastic(rf, props)

# the model -> rh rf = VRHAvg(["VSAND", "VCLAY"], [sand, clay]) # check modulus variation with volume fraction props = dict(VCLAY=np.linspace(0, 1, 51),) plot_elastic(rf, props)

Apply Nur's critical porosity adjustment to turn average bulk modulus in a porous frame modulus.

In [5]:

pn = NurCriticalPoroAdjust(["poro"], rf, 0.39)
props = dict(
    VCLAY=np.linspace(0, 1, 51),
    poro=0.2,
)
plot_elastic(pn, props)

pn = NurCriticalPoroAdjust(["poro"], rf, 0.39) props = dict( VCLAY=np.linspace(0, 1, 51), poro=0.2, ) plot_elastic(pn, props)

Create a StressModel and apply MacBeths Stress Adjustment

In [6]:

grad = 1/145.038*3.28084 # 1 psi/ft converted to MPa/m
# reference stress to surface
stressm = LGStressModel(grad, 0, 0)
# Then apply the MacBeth adjustment model
sm = MacBethStressAdjust(["depth, pres_init"], pn, stressm, 0.2, 17, 0.4, 18)
props = dict(
    VCLAY=np.linspace(0, 1, 51),
    poro=0.2,
    pres_init=60, #MPa
    depth = 3500, #m
    pres=40 # 20MPa pressure decline
)
plot_elastic(sm, props)

grad = 1/145.038*3.28084 # 1 psi/ft converted to MPa/m # reference stress to surface stressm = LGStressModel(grad, 0, 0) # Then apply the MacBeth adjustment model sm = MacBethStressAdjust(["depth, pres_init"], pn, stressm, 0.2, 17, 0.4, 18) props = dict( VCLAY=np.linspace(0, 1, 51), poro=0.2, pres_init=60, #MPa depth = 3500, #m pres=40 # 20MPa pressure decline ) plot_elastic(sm, props)

In [7]:

props = dict(
    VCLAY=np.linspace(0, 1, 2),
    poro=0.2,
    pres_init=60, #MPa
    depth = 3500, #m
    pres=40 # 20MPa depressure
)
sm.trace_tree(props, ["bulk_modulus", "vp"])

props = dict( VCLAY=np.linspace(0, 1, 2), poro=0.2, pres_init=60, #MPa depth = 3500, #m pres=40 # 20MPa depressure ) sm.trace_tree(props, ["bulk_modulus", "vp"])

Out[7]:

MacBethStressAdjust_140377422097856
├── bulk_modulus : [15.20625479 11.18854867]
├── vp : [3795.38603058 2980.2424466 ]
└── NurCriticalPoroAdjust_140378536950080
    ├── bulk_modulus : [14.56666667 10.71794872]
    ├── vp : [3675.62361061 2894.5555797 ]
    └── VRHAvg_140378536698880
        ├── bulk_modulus : [29.9 22. ]
        ├── vp : [4710.11519872 3709.21826438]
        ├── Mineral_140378536987760
        │   ├── bulk_modulus : 29.9
        │   └── vp : 4710.1151987170315
        └── Mineral_140378536988432
            ├── bulk_modulus : 22
            └── vp : 3709.2182643789142

Create a Hashin-Shtrikman-Walpole Model¶

In [8]:

import numpy as np
from digirock import RockFrame, Element, Mineral, WaterBW92, HSWAvg
    
sand = Mineral(2.7, 35, 45)
clay = Mineral(2.51, 75, 31)
water = WaterBW92()
# water.vp = water.velocity
# water.vs = lambda x: 0
hs = HSWAvg(["VSAND", "VCLAY", "poro"], [sand, clay, water])

# use xarray to create an orthogonal set of properties in VCLAY and POROSITY
ds, props = create_orthogonal_props(VCLAY=np.linspace(0.1, 0.25, 16), poro=np.linspace(0.05, 0.3, 20), temp=np.r_[20], pres=np.r_[10])
ds["bulk_modulus"] = (ds.dims, hs.bulk_modulus(props))
ds["shear_modulus"] = (ds.dims, hs.shear_modulus(props))

fig, axs = plt.subplots(ncols=2, figsize=figsize)
ds.sel(temp=20, pres=10).bulk_modulus.plot(ax=axs[0])
ds.sel(temp=20, pres=10).shear_modulus.plot(ax=axs[1])

import numpy as np from digirock import RockFrame, Element, Mineral, WaterBW92, HSWAvg sand = Mineral(2.7, 35, 45) clay = Mineral(2.51, 75, 31) water = WaterBW92() # water.vp = water.velocity # water.vs = lambda x: 0 hs = HSWAvg(["VSAND", "VCLAY", "poro"], [sand, clay, water]) # use xarray to create an orthogonal set of properties in VCLAY and POROSITY ds, props = create_orthogonal_props(VCLAY=np.linspace(0.1, 0.25, 16), poro=np.linspace(0.05, 0.3, 20), temp=np.r_[20], pres=np.r_[10]) ds["bulk_modulus"] = (ds.dims, hs.bulk_modulus(props)) ds["shear_modulus"] = (ds.dims, hs.shear_modulus(props)) fig, axs = plt.subplots(ncols=2, figsize=figsize) ds.sel(temp=20, pres=10).bulk_modulus.plot(ax=axs[0]) ds.sel(temp=20, pres=10).shear_modulus.plot(ax=axs[1])

Out[8]:

<matplotlib.collections.QuadMesh at 0x7fac29ae3ee0>

Create a Cemented Sand Model¶

In [9]:

from digirock import CementedSand, Mineral
import numpy as np
import xarray as xr

sand = Mineral(2.7, 35, 45)
cement = Mineral(2.9, 55, 50)

cs = CementedSand("VSAND", sand, "VCEMENT", cement, ncontacts=9)
cs.bulk_modulus({"poro":0.2, "VCEMENT": 0.05, "ncontacts":40})

from digirock import CementedSand, Mineral import numpy as np import xarray as xr sand = Mineral(2.7, 35, 45) cement = Mineral(2.9, 55, 50) cs = CementedSand("VSAND", sand, "VCEMENT", cement, ncontacts=9) cs.bulk_modulus({"poro":0.2, "VCEMENT": 0.05, "ncontacts":40})

Out[9]:

19.88948478318842

In [10]:

ds

ds

Out[10]:

<xarray.Dataset>
Dimensions:        (VCLAY: 16, poro: 20, temp: 1, pres: 1)
Coordinates:
  * VCLAY          (VCLAY) float64 0.1 0.11 0.12 0.13 ... 0.22 0.23 0.24 0.25
  * poro           (poro) float64 0.05 0.06316 0.07632 ... 0.2737 0.2868 0.3
  * temp           (temp) int64 20
  * pres           (pres) int64 10
Data variables:
    bulk_modulus   (VCLAY, poro, temp, pres) float64 28.1 26.7 ... 17.57 17.16
    shear_modulus  (VCLAY, poro, temp, pres) float64 19.66 19.16 ... 11.36 11.03

xarray.Dataset

• Dimensions:
  • VCLAY: 16
  • poro: 20
  • temp: 1
  • pres: 1
• Coordinates: (4)
  • VCLAY
    (VCLAY)
    float64
    0.1 0.11 0.12 ... 0.23 0.24 0.25

    array([0.1 , 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2 , 0.21,
           0.22, 0.23, 0.24, 0.25])

  • poro
    (poro)
    float64
    0.05 0.06316 0.07632 ... 0.2868 0.3

    array([0.05    , 0.063158, 0.076316, 0.089474, 0.102632, 0.115789, 0.128947,
           0.142105, 0.155263, 0.168421, 0.181579, 0.194737, 0.207895, 0.221053,
           0.234211, 0.247368, 0.260526, 0.273684, 0.286842, 0.3     ])

  • temp
    (temp)
    int64
    20

    array([20])

  • pres
    (pres)
    int64
    10

    array([10])

• Data variables: (2)
  • bulk_modulus
    (VCLAY, poro, temp, pres)
    float64
    28.1 26.7 25.5 ... 17.57 17.16

    array([[[[28.09716567]],
    
            [[26.69544647]],
    
            [[25.4982223 ]],
    
            [[24.4542517 ]],
    
            [[23.52836821]],
    
            [[22.69562705]],
    
            [[21.93783616]],
    
            [[21.24141033]],
    
            [[20.59599495]],
    
            [[19.9935555 ]],
    
    ...
    
            [[21.47733675]],
    
            [[20.90658521]],
    
            [[20.3653904 ]],
    
            [[19.85009797]],
    
            [[19.35769591]],
    
            [[18.8856772 ]],
    
            [[18.43193632]],
    
            [[17.99469028]],
    
            [[17.57241765]],
    
            [[17.16381096]]]])

  • shear_modulus
    (VCLAY, poro, temp, pres)
    float64
    19.66 19.16 18.66 ... 11.36 11.03

    array([[[[19.6601088 ]],
    
            [[19.15610639]],
    
            [[18.66373305]],
    
            [[18.18259088]],
    
            [[17.71229992]],
    
            [[17.25249719]],
    
            [[16.8028357 ]],
    
            [[16.36298359]],
    
            [[15.93262334]],
    
            [[15.511451  ]],
    
    ...
    
            [[14.24633227]],
    
            [[13.86021559]],
    
            [[13.48191364]],
    
            [[13.11119154]],
    
            [[12.74782375]],
    
            [[12.39159357]],
    
            [[12.04229274]],
    
            [[11.69972102]],
    
            [[11.36368583]],
    
            [[11.03400185]]]])

• Attributes: (0)

In [11]:

ds, props = create_orthogonal_props(
    VCEMENT=np.linspace(0, 0.15, 16), ncontacts=np.arange(10, 100, 20), poro=np.arange(0.05, 0.31, 0.01)
)

ds["bulk_modulus"] = (ds.dims, cs.bulk_modulus(props))
fig, axs = plt.subplots(ncols=3, figsize=figsize)
ds.sel(poro=0.1).bulk_modulus.plot(ax=axs[0], vmin=0, vmax=50)
ds.sel(poro=0.2).bulk_modulus.plot(ax=axs[1], vmin=0, vmax=50)
# or slice in the porosity direction
ds.sel(ncontacts=50).bulk_modulus.plot(ax=axs[2], vmin=0, vmax=50)
fig.tight_layout()

ds, props = create_orthogonal_props( VCEMENT=np.linspace(0, 0.15, 16), ncontacts=np.arange(10, 100, 20), poro=np.arange(0.05, 0.31, 0.01) ) ds["bulk_modulus"] = (ds.dims, cs.bulk_modulus(props)) fig, axs = plt.subplots(ncols=3, figsize=figsize) ds.sel(poro=0.1).bulk_modulus.plot(ax=axs[0], vmin=0, vmax=50) ds.sel(poro=0.2).bulk_modulus.plot(ax=axs[1], vmin=0, vmax=50) # or slice in the porosity direction ds.sel(ncontacts=50).bulk_modulus.plot(ax=axs[2], vmin=0, vmax=50) fig.tight_layout()

In [ ]: