<jupyter_start><jupyter_text>Gassmann Fluid SubstitutionThis notebook is based [on the original by <NAME>](https://github.com/ezygeo-ai/Fluid-Substitution), which is licenced under the MIT licence.Original: Copyright (C) 2020 ezygeo.com, <NAME> (<EMAIL>)Modified: <NAME>, 2021. Same licence. I am going to check the result using `bruges`.<jupyter_code>import bruges as bg<jupyter_output><empty_output><jupyter_text>This Notebook requires bruges 0.4.0 or greater.<jupyter_code>bg.__version__
import numpy as np<jupyter_output><empty_output><jupyter_text>DataThe original notebook had cgs units. I'm going to do everything in SI because science.<jupyter_code># Converting to SI units.
Vp = 3.349 * 1000
Vs = 1.836 * 1000

So = 0.7  # Initial condition.
Sw = 1    # Subbed case.

Por = 0.3

Fq = 0.7
Fc = 0.3
Kq = 36 * 1e9
Kc = 75 * 1e9
Ko = 1.6 * 1e9
Kw = 2.83 * 1e9

Rho = 2.13 * 1000
Rhoo = 0.8 * 1000
Rhow = 1 * 1000<jupyter_output><empty_output><jupyter_text>Define FormulaeThe ezygeo functions all used rounding on the outputs... I have removed this 'feature'.  1. K Saturated Initial _Biot (1956)_<jupyter_code>def msat(vs, rho):
    '''
    Miu saturated (C) hpratama 2020
    
    Input:
    vs = Shear Velocity (Km/s)
    rho = Density (g/cc)
  
    '''
    return rho*(vs**2)<jupyter_output><empty_output><jupyter_text>I'm going to check all these formula against the equivalent in `bruges`, if it's in there.<jupyter_code>assert msat(Vs, Rho) == bg.rockphysics.moduli.mu(vs=Vs, rho=Rho)
def ksat1(vp, vs, rho):
    '''
    K saturated (C) hpratama 2020
    
    Input:
    vp = P wave velocity (Km/s)
    vs = Shear wave Velocity (Km/s)
    rho = Density (g/cc)
  
    '''
    return ((vp**2)*rho) - 4/3*(msat(vs, rho))
assert ksat1(Vp, Vs, Rho) == bg.rockphysics.moduli.bulk(Vp, Vs, Rho)<jupyter_output><empty_output><jupyter_text>2. K Fluid  _Batzle and Wang, 1992_<jupyter_code>def kfl(sw, khc, kw):
    '''
    K Fluid (C) hpratama 2020
    
    Input:
    sw = Water Saturation
    khc = K Hydrocarbon
    kw = K Water
    
    '''
    return ((sw/kw) + ((1-sw)/khc))**-1
ezygeo = kfl(1 - So, Ko, Kw)
bruges = bg.rockphysics.fluids.wood(Kw, Ko, 1 - So)

assert ezygeo == bruges<jupyter_output><empty_output><jupyter_text>3. K Mineral _Voight-Reuss-Hill Average_<jupyter_code>def km(por, fq, fc, kq, kc):
    '''
    K mineral (K0) (C) hpratama 2020
    
    Input:
    por = Porosity
    fq = Quartz Fraction
    fc = Calcite Fraction
    kq = Quarzt Bulk
    kc = Calcite Bulk
    '''
    kv = ((fq)*kq) + ((fc)*kc) 
    kr = (((fq)/kq) + ((fc)/kc))**(-1)
    return (kv + kr)/2
ezygeo = km(Por, Fq, Fc, Kq, Kc)
bruges = bg.rockphysics.vrh(Kc, Kq, Fc)

assert ezygeo == bruges<jupyter_output><empty_output><jupyter_text>4. K Saturated Final  _Mavko et al._<jupyter_code>def ksat2(ksat1, k0, kfl1, kfl2, por):
    '''
    K Saturation after fluid substitution (C) hpratama 2020
    
    Input:
    ksat1 = K Saturation before Fluid Substitution
    k0 = K Mineral
    kfl1 = K Fluid before Substitution
    kfl2 = K Fluid after Substitution
    por = Porosity
    
    '''
    a = (ksat1/(k0 - ksat1)) - (kfl1/(por*(k0 - kfl1))) + (kfl2/(por*(k0 - kfl2)))
    return (a*(k0))/(1+a)
# First Ezygeo's solution.
Ksat1 = ksat1(Vp, Vs, Rho)
Kfl1 = kfl((1-So), Ko, Kw)
Kfl2 = kfl(Sw, Ko, Kw)
K0 = km(Por, Fq, Fc, Kq, Kc)
Ksat2_ezygeo = ksat2(Ksat1, K0, Kfl1, Kfl2, Por)

# And using bruges.
Ksat1 = bg.rockphysics.moduli.bulk(Vp, Vs, Rho)
Kfl1 = bg.rockphysics.fluids.wood(Kw, Ko, 1 - So)
Kfl2 = bg.rockphysics.fluids.wood(Kw, Ko, Sw)
K0 = bg.rockphysics.vrh(Kc, Kq, Fc)
Ksat2_bruges = bg.rockphysics.fluidsub.avseth_gassmann(Ksat1, Kfl1, Kfl2, K0, Por)

np.testing.assert_almost_equal(Ksat2_ezygeo, Ksat2_bruges, decimal=4)<jupyter_output><empty_output><jupyter_text>5. Vp   Compute fluid density:<jupyter_code>def rhofluid(sw, so, rhoo, rhow):
    
    rho1 = ((1-so) * rhow) + (so * rhoo)
    rho2 = ((sw) * rhow) + ((1 - sw) * rhoo)
    return rho1, rho2
# Lines adapted from bruges source code.
rhofl = (1 - So) * Rhow + So * Rhoo
rhofl2 = 1 * Rhow

assert rhofluid(Sw, So, Rhoo, Rhow) == (rhofl, rhofl2)<jupyter_output><empty_output><jupyter_text>Compute new rock density:<jupyter_code>def rhosat(rho, por, rhow, sw, so, rhoo):
    rho1, rho2 = rhofluid(sw, so, rhoo, rhow)
    return rho + (por * (rho2 - rho1))
# Lines adapted from bruges source code.
rhog = (Rho - Por * rhofl) / (1 - Por)
rho2 = Por * rhofl2 + (1 - Por) * rhog

# ezygeo solution.
Rhosat = rhosat(Rho, Por, Rhow, Sw, So, Rhoo)

assert Rhosat == rho2<jupyter_output><empty_output><jupyter_text>Compute new Vp:<jupyter_code>def vp(ksat2, vs, rhosat):
    return ((ksat2 + (4/3*(msat(vs, rhosat)))) / (rhosat))**0.5
ezygeo = vp(Ksat2_ezygeo, Vs, Rhosat)
#                         ^^ Error: Vs(2) != Vs(1)

G = bg.rockphysics.moduli.mu(vs=Vs, rho=Rho)
bruges = bg.rockphysics.moduli.vp(bulk=Ksat2_bruges, mu=G, rho=rho2)

ezygeo, bruges<jupyter_output><empty_output><jupyter_text>Calculate `ezygeo` result To get new value of Vp, there is some variable that need to calculate.* First we need to calculate K saturation before fluid subtitution. This K saturartion could be calculated from the old<jupyter_code>Ksat1 = ksat1(Vp, Vs, Rho)<jupyter_output><empty_output><jupyter_text>* Then, we estimate the K fluid before and after fluid substitution<jupyter_code>Kfl1 = kfl((1-So), Ko, Kw)
Kfl2 = kfl(Sw, Ko, Kw)  # = Kw<jupyter_output><empty_output><jupyter_text>* After that, we evaluate K Mineral that can be estimate from Hill formula, this formula need Voight-Reuss Bounds calculation too. First we need to calculate K saturation before fluid subtitution.<jupyter_code>K0 = km(Por, Fq, Fc, Kq, Kc)<jupyter_output><empty_output><jupyter_text>* Here we calculate K saturated after fluid substitution<jupyter_code>Ksat2 = ksat2(Ksat1, K0, Kfl1, Kfl2, Por)<jupyter_output><empty_output><jupyter_text>* To calculate the new Vp, we need to calculate density saturation. After that we estimate new Vp with new K saturated after fluid substitution<jupyter_code>Rhosat = rhosat(Rho, Por, Rhow, Sw, So, Rhoo)
Vp2 = vp(Ksat2, Vs, Rhosat)
print(Vp2)<jupyter_output>3439.9002222704626<jupyter_text>Done.  Compute using Avseth's algorithm in `bruges`<jupyter_code>from bruges.rockphysics.fluidsub import smith_gassmann, avseth_gassmann, avseth_fluidsub, smith_fluidsub
avseth_fluidsub(Vp, Vs, Rho,
                Por,
                rhof1=(1 - So) * Rhow + So * Rhoo,
                rhof2=Rhow,
                kmin=bg.rockphysics.vrh(Kc, Kq, Fc),
                kf1=bg.rockphysics.fluids.wood(Kw, Ko, 1-So),
                kf2=Kw,
               )<jupyter_output><empty_output><jupyter_text>Compute using Smith's algorithm in `bruges`<jupyter_code>smith_fluidsub(Vp, Vs, Rho,
               Por,
               Rhow, Rhoo,
               sw=1 - So, swnew=1,  # Fluid 1: Mix, Fluid 2: Brine
               kw=Kw,
               khc=Ko,
               kclay=Kc,
               kqtz=Kq,
               vclay=Fc,
               )<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13072546230440968, 'fraction_numerical': 0.033285917496443815, 'mean_word_length': 3.539057456423499, 'pattern_counts': {'":': 0, '<': 76, '<?xml version=': 0, '>': 76, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '616733', 'license': 'Apache-2.0', 'n_tokens_mistral': 3080, 'n_tokens_neox': 2881, 'n_words': 712, 'path': 'notebooks_dev/Fluidsubs_from_ezygeo-ai_for_checking.ipynb', 'repo_name': 'agile-geoscience/in-bruges', 'size': '7071'}