Skip to content

Models

Methods can be imported from digirock.models.

Basic Models

Functions related to the moduli and mixed moduli of minerals

This module contains functions related to the calculation of

  1. modulus bounds
  2. moduli for mixed mineral materials
  3. Gassmann fluid substitution

Refs

MacBeth, 2004, A classification for the pressure-sensitivity properties of a sandstone rock frame

Smith, Sondergeld and Rai, 2003, Gassmann fluid substitutions: A tutorial, Geophysics, 68, pp430-440

dryframe_acoustic(ksat, kfl, k0, phi)

Dry frame bulk modulus \(K^*\) from material saturated modulus

This is essentially inverse Gassmann fluid substitution and assumes you know the bulk modulus of the matrix material k0.

\[ K^* = \frac{K_{sat}\left( \frac{\phi K_0}{K_{fl}} + 1 - \phi \right) - K_0}{\frac{\phi K_0}{K_{fl}}+\frac{K_{sat}}{K_0}-1-\phi} \]

Parameters:

Name Type Description Default
ksat Union[numpy.ndarray, float]

saturated bulk modulus

 required
kfl Union[numpy.ndarray, float]

fluid bulk modulus

 required
k0 Union[numpy.ndarray, float]

matrix bulk modulus

 required
phi Union[numpy.ndarray, float]

porosity (frac)

 required

Returns:

Type Description
Union[numpy.ndarray, float]

backed-out dry fame bulk modulus
Source code in digirock/models/_mod.py
Python
def dryframe_acoustic(
    ksat: NDArrayOrFloat, kfl: NDArrayOrFloat, k0: NDArrayOrFloat, phi: NDArrayOrFloat
) -> NDArrayOrFloat:
    """Dry frame bulk modulus $K^*$ from material saturated modulus

    This is essentially inverse Gassmann fluid substitution and assumes you know the bulk modulus of
    the matrix material `k0`.

    $$
    K^* = \\frac{K_{sat}\\left( \\frac{\\phi K_0}{K_{fl}} + 1 - \\phi \\right) - K_0}{\\frac{\\phi K_0}{K_{fl}}+\\frac{K_{sat}}{K_0}-1-\\phi}
    $$

    Args:
        ksat: saturated bulk modulus
        kfl: fluid bulk modulus
        k0: matrix bulk modulus
        phi: porosity (frac)

    Returns:
        backed-out dry fame bulk modulus
    """
    return (ksat * (phi * k0 / kfl + 1 - phi) - k0) / (
        phi * k0 / kfl + ksat / k0 - 1 - phi
    )

dryframe_density(rhob, rhofl, phi)

Dry frame grain density

Primarily used for fluid substitution.

Parameters:

Name Type Description Default
rhob array-like

saturated rock bulk modulus

 required
rhofl array-like

fluid density

 required
phi array-like

porosity (frac)

 required

Returns:

Type Description
(array-like)

dry frame grain density
Source code in digirock/models/_mod.py
Python
def dryframe_density(rhob, rhofl, phi):
    """Dry frame grain density

    Primarily used for fluid substitution.

    Args:
        rhob (array-like): saturated rock bulk modulus
        rhofl (array-like): fluid density
        phi (array-like): porosity (frac)

    Returns:
        (array-like): dry frame grain density
    """
    return (rhob - rhofl * phi) / (1 - phi)

dryframe_dpres(dry_mod, pres1, pres2, sse, ssp)

Calculates the dry rock frame for a given stress regime by factoring the difference between two formation pressure scenarios.

\[ m(P) = m \frac{1 + S_E e^{\tfrac{-P_1}{S_P}}}{1 + S_E e^{\tfrac{-P_2}{S_P}}} \]

Parameters:

Name Type Description Default
dry_mod Union[numpy.ndarray, float]

Dryframe modulus (MPa)

 required
pres1 Union[numpy.ndarray, float]

Pressure calibrated to dryframe modulus. (MPa)

 required
pres2 Union[numpy.ndarray, float]

Pressure of output (MPa)

 required
sse Union[numpy.ndarray, float]

Stress-sensitivity parameter E

 required
ssp Union[numpy.ndarray, float]

Stress-sensitivity parameter P

 required

Returns:

Type Description
Union[numpy.ndarray, float]

the dry-frame modulus adjusted to pressure (pres2)

References

[1] Amini and Alvarez (2014) [2] MacBeth (2004)

Source code in digirock/models/_mod.py
Python
def dryframe_dpres(
    dry_mod: NDArrayOrFloat,
    pres1: NDArrayOrFloat,
    pres2: NDArrayOrFloat,
    sse: NDArrayOrFloat,
    ssp: NDArrayOrFloat,
) -> NDArrayOrFloat:
    """Calculates the dry rock frame for a given stress regime by factoring the
    difference between two formation pressure scenarios.

    $$
    m(P) = m \\frac{1 + S_E e^{\\tfrac{-P_1}{S_P}}}{1 + S_E e^{\\tfrac{-P_2}{S_P}}}
    $$

    Args:
        dry_mod: Dryframe modulus (MPa)
        pres1: Pressure calibrated to dryframe modulus. (MPa)
        pres2: Pressure of output (MPa)
        sse: Stress-sensitivity parameter E
        ssp: Stress-sensitivity parameter P

    Returns:
        the dry-frame modulus adjusted to pressure (pres2)

    References:
        [1] Amini and Alvarez (2014)
        [2] MacBeth (2004)
    """
    return (
        dry_mod
        * (1 + (sse * np.exp(-pres1 / ssp)))
        / (1 + (sse * np.exp(-pres2 / ssp)))
    )

dryframe_stress(mod_e, mod_p, inf, p)

Dry frame bulk modulus from stress sensitivity coefficients.

Based on MacBeth 2004 for sandstone rock frame.

Parameters:

Name Type Description Default
mod_e array-like

modulus stress sensitivity metric *2 MPa

 required
mod_p array-like

modulus characteristic pressure constant *2 MPa

 required
inf array-like

modulus background high pressure asymptote MPa

 required
p array-like

the effective pressure MPa

 required
Source code in digirock/models/_mod.py
Python
def dryframe_stress(mod_e, mod_p, inf, p):
    """Dry frame bulk modulus from stress sensitivity coefficients.

    Based on MacBeth 2004 for sandstone rock frame.

    Args:
        mod_e (array-like): modulus stress sensitivity metric *2 MPa
        mod_p (array-like): modulus characteristic pressure constant *2 MPa
        inf (array-like): modulus background high pressure asymptote MPa
        p (array-like): the effective pressure MPa
    """
    return inf / (1 + (mod_e * np.exp(-p / mod_p)))

gassmann_fluidsub(kdry, kfl, k0, phi)

Gassmann fluid substitution for saturated rock bulk modulus

Gassmann fluid substitution assumes: The rock is homogenous and isotropic, rocks with diverse mineral acoustic properties will violate this assumption. All the porespace is connected, this is propably not true for low porosity rocks. The frequencies being used to measure the rock properties fall within the seismic bandwidth. Logging frequencies are usually ok for high-porosity clean sands. Homogeonous mixing of the fluid. This is reasonable over geologic time but may be inappropriate for production timescales.

Parameters:

Name Type Description Default
kdry Union[numpy.ndarray, float]

dry frame rock bulk modulus

 required
kfl Union[numpy.ndarray, float]

fluid bulk modulus

 required
k0 Union[numpy.ndarray, float]

matrix bulk modulus

 required
phi Union[numpy.ndarray, float]

porosity (frac)

 required

Returns:

Type Description
Union[numpy.ndarray, float]

ksat saturated rock bulk modulus
Source code in digirock/models/_mod.py
Python
def gassmann_fluidsub(
    kdry: NDArrayOrFloat, kfl: NDArrayOrFloat, k0: NDArrayOrFloat, phi: NDArrayOrFloat
) -> NDArrayOrFloat:
    """Gassmann fluid substitution for saturated rock bulk modulus

    Gassmann fluid substitution assumes:
        The rock is homogenous and isotropic, rocks with diverse mineral acoustic properties will
            violate this assumption.
        All the porespace is connected, this is propably not true for low porosity rocks.
        The frequencies being used to measure the rock properties fall within the seismic
            bandwidth. Logging frequencies are usually ok for high-porosity clean sands.
        Homogeonous mixing of the fluid. This is reasonable over geologic time but may be
            inappropriate for production timescales.

    Args:
        kdry: dry frame rock bulk modulus
        kfl: fluid bulk modulus
        k0: matrix bulk modulus
        phi: porosity (frac)

    Returns:
        ksat saturated rock bulk modulus
    """
    return kdry + np.power(1 - kdry / k0, 2) / (
        safe_divide(phi, kfl)
        + safe_divide((1 - phi), k0)
        - safe_divide(kdry, np.power(k0, 2))
    )

reuss_lower_bound(mod, frac, *argv)

Calculate the Reuss lower bound for material moduli

Parameters:

Name Type Description Default
mod Union[numpy.ndarray, float]

first fluid modulus

 required
frac Union[numpy.ndarray, float]

first fluid volume fraction

 required
argv Union[numpy.ndarray, float]

additional fluid and volume fraction pairs.

 ()

Returns:

Type Description
Union[numpy.ndarray, float]

Reuss upper bound for modulus mix
Source code in digirock/models/_mod.py
Python
def reuss_lower_bound(
    mod: NDArrayOrFloat, frac: NDArrayOrFloat, *argv: NDArrayOrFloat
) -> NDArrayOrFloat:
    """Calculate the Reuss lower bound for material moduli

    Args:
        mod: first fluid modulus
        frac: first fluid volume fraction
        argv: additional fluid and volume fraction pairs.

    Returns:
        Reuss upper bound for modulus mix
    """
    args = _process_vfrac(*((mod, frac) + argv))
    mod_sum = 0
    for modi, vfraci in chunked(args, 2):
        mod_sum = mod_sum + safe_divide(np.array(vfraci), np.array(modi))
    return safe_divide(1, mod_sum)

saturated_density(rhog, rhofl, phi)

Saturated bulk density

Parameters:

Name Type Description Default
rhog array-like

dry frame grain density

 required
rhofl array-like

fluid density

 required
phi array-like

porosity (frac)

 required

Returns:

Type Description
(array-like)

saturated rock bulk density (rhob)
Source code in digirock/models/_mod.py
Python
def saturated_density(rhog, rhofl, phi):
    """Saturated bulk density

    Args:
        rhog (array-like): dry frame grain density
        rhofl (array-like): fluid density
        phi (array-like): porosity (frac)

    Returns:
        (array-like): saturated rock bulk density (rhob)
    """
    return rhog * (1 - phi) + rhofl * phi

voigt_upper_bound(mod, frac, *argv)

Calculate the Voigt upper bound for material moduli

Parameters:

Name Type Description Default
mod Union[numpy.ndarray, float]

first fluid modulus

 required
frac Union[numpy.ndarray, float]

first fluid volume fraction

 required
argv Union[numpy.ndarray, float]

additional fluid and volume fraction pairs.

 ()

Returns:

Type Description
Union[numpy.ndarray, float]

Voigt upper bound for modulus mix
Source code in digirock/models/_mod.py
Python
def voigt_upper_bound(
    mod: NDArrayOrFloat, frac: NDArrayOrFloat, *argv: NDArrayOrFloat
) -> NDArrayOrFloat:
    """Calculate the Voigt upper bound for material moduli

    Args:
        mod: first fluid modulus
        frac: first fluid volume fraction
        argv: additional fluid and volume fraction pairs.

    Returns:
        Voigt upper bound for modulus mix
    """
    args = _process_vfrac(*((mod, frac) + argv))
    mod_sum = 0
    for modi, vfraci in chunked(args, 2):
        mod_sum = mod_sum + np.array(vfraci) * np.array(modi)
    return mod_sum

vrh_avg(mod, frac, *argv)

Calculates the Voigt-Reuss-Hill average for a multi-modulus mix.

Parameters:

Name Type Description Default
mod Union[numpy.ndarray, float]

first modulus

 required
frac Union[numpy.ndarray, float]

first volume fraction

 required
argv Union[numpy.ndarray, float]

additional modulus and volume fraction pairs.

 ()

Returns:

Type Description
Union[numpy.ndarray, float]

Voigt-Reuss-Hill average modulus
Source code in digirock/models/_mod.py
Python
def vrh_avg(
    mod: NDArrayOrFloat, frac: NDArrayOrFloat, *argv: NDArrayOrFloat
) -> NDArrayOrFloat:
    """Calculates the Voigt-Reuss-Hill average for a multi-modulus mix.

    Args:
        mod: first modulus
        frac: first volume fraction
        argv: additional modulus and volume fraction pairs.

    Returns:
        Voigt-Reuss-Hill average modulus
    """
    vub = voigt_upper_bound(mod, frac, *argv)
    rlb = reuss_lower_bound(mod, frac, *argv)
    return 0.5 * (vub + rlb)

Hashin-Shtrikman Model

Hashin-Shtrikman-Walpole Models for moduli mixing.

Refs

Mavko et. al. (2009) The Rock Physics Handbook, Cambridge Press

hs_kbounds2(k1, k2, mu1, f1)

Calculate the Hashin-Shtrikman bulk moduli bounds for an isotropic elastic mixture

Upper and lower bounds are computed by interchanging which constituent material is subscripted 1 and which is subscripted 2. In general:

  • Upper bound -> stiffest material is subscripted 1
  • Lower bound -> stiffest material is subscripted 1

Parameters:

Name Type Description Default
k1 Union[numpy.ndarray, float]

constituent one bulk modulus

 required
k2 Union[numpy.ndarray, float]

constituent two bulk modulus

 required
mu1 Union[numpy.ndarray, float]

constituent one shear modulus

 required
f1

material one volume fraction 0 <= f1 <= 1

 required

Returns:

Type Description
Union[numpy.ndarray, float]

HS bulk modulus bounds

Refs

Rock Physics Handbook p170

Source code in digirock/models/_hsw.py
Python
def hs_kbounds2(
    k1: NDArrayOrFloat, k2: NDArrayOrFloat, mu1: NDArrayOrFloat, f1
) -> NDArrayOrFloat:
    """Calculate the Hashin-Shtrikman bulk moduli bounds for an isotropic elastic mixture

    Upper and lower bounds are computed by interchanging which constituent material is
    subscripted 1 and which is subscripted 2. In general:

     - Upper bound -> stiffest material is subscripted 1
     - Lower bound -> stiffest material is subscripted 1

    Args:
        k1: constituent one bulk modulus
        k2: constituent two bulk modulus
        mu1: constituent one shear modulus
        f1: material one volume fraction 0 <= f1 <= 1

    Returns:
        HS bulk modulus bounds

    Refs:
        Rock Physics Handbook p170
    """
    return k1 + (1 - f1) / (1 / (k2 - k1) + f1 / (k1 + 4 * mu1 / 3))

hs_mubounds2(k1, mu1, mu2, f1)

Calculate the Hashin-Shtrikman shear moduli bounds for an isotropic elastic mixture

Upper and lower bounds are computed by interchanging which constituent material is subscripted 1 and which is subscripted 2. In general:

  • Upper bound -> stiffest material is subscripted 1
  • Lower bound -> stiffest material is subscripted 1

Parameters:

Name Type Description Default
k1 Union[numpy.ndarray, float]

constituent one bulk modulus

 required
mu1 Union[numpy.ndarray, float]

constituent one shear modulus

 required
mu2 Union[numpy.ndarray, float]

constituent two shear modulus

 required
f1 Union[numpy.ndarray, float]

material one volume fraction 0 <= f1 <= 1

 required

Returns:

Type Description
Union[numpy.ndarray, float]

HS shear modulus bounds (positive)

Refs

Rock Physics Handbook p170

Source code in digirock/models/_hsw.py
Python
def hs_mubounds2(
    k1: NDArrayOrFloat, mu1: NDArrayOrFloat, mu2: NDArrayOrFloat, f1: NDArrayOrFloat
) -> NDArrayOrFloat:
    """Calculate the Hashin-Shtrikman shear moduli bounds for an isotropic elastic mixture

    Upper and lower bounds are computed by interchanging which constituent material is
    subscripted 1 and which is subscripted 2. In general:

     - Upper bound -> stiffest material is subscripted 1
     - Lower bound -> stiffest material is subscripted 1

    Args:
        k1: constituent one bulk modulus
        mu1: constituent one shear modulus
        mu2: constituent two shear modulus
        f1: material one volume fraction 0 <= f1 <= 1

    Return:
        HS shear modulus bounds (positive)

    Refs:
        Rock Physics Handbook p170
    """
    return mu1 + (1 - f1) / (
        (1 / (mu2 - mu1)) + 2 * f1 * (k1 + 2 * mu1) / (5 * mu1 * (k1 + 4 * mu1 / 3))
    )

hsw_avg(bulk_mod, shear_mod, frac, *argv)

Hashin-Shtrikman-Walpole average moduli for N materials.

Parameters:

Name Type Description Default
bulk_mod Union[numpy.ndarray, float]

first component bulk modulus

 required
shear_mod Union[numpy.ndarray, float]

first component shear modulus

 required
frac Union[numpy.ndarray, float]

first component volume fraction

 required
*argv Union[numpy.ndarray, float]

additional triplets of bulk_mod, shear_mod and frac for additional components

 ()

Returns:

Type Description
Tuple[Union[numpy.ndarray, float], Union[numpy.ndarray, float]]

K_hs_avg, Mu_hs_avg
Source code in digirock/models/_hsw.py
Python
def hsw_avg(
    bulk_mod: NDArrayOrFloat,
    shear_mod: NDArrayOrFloat,
    frac: NDArrayOrFloat,
    *argv: NDArrayOrFloat,
) -> Tuple[NDArrayOrFloat, NDArrayOrFloat]:
    """Hashin-Shtrikman-Walpole average moduli for N materials.

    Args:
        bulk_mod: first component bulk modulus
        shear_mod: first component shear modulus
        frac: first component volume fraction
        *argv: additional triplets of bulk_mod, shear_mod and frac for additional components

    Returns:
        K_hs_avg, Mu_hs_avg
    """
    k_hsp, k_hsm, mu_hsp, mu_hsm = hsw_bounds(bulk_mod, shear_mod, frac, *argv)
    return 0.5 * (k_hsp + k_hsm), 0.5 * (mu_hsp + mu_hsm)

hsw_bounds(bulk_mod, shear_mod, frac, *argv)

Hashin-Shtrikman-Walpole bounds

\[ K^{HS^+} = \Lambda(\mu_{max}) \]

Parameters:

Name Type Description Default
bulk_mod Union[numpy.ndarray, float]

first component bulk modulus

 required
shear_mod Union[numpy.ndarray, float]

first component shear modulus

 required
frac Union[numpy.ndarray, float]

first component volume fraction

 required
*argv Union[numpy.ndarray, float]

additional triplets of bulk_mod, shear_mod and frac for additional components

 ()

Returns:

Type Description
Tuple[Union[numpy.ndarray, float], Union[numpy.ndarray, float], Union[numpy.ndarray, float], Union[numpy.ndarray, float]]

K_hsp, K_hsm, Mu_hsp, Mu_hsm

Refs

Rock Physics Handbook p171-172

This example comes form the RPH p172

Examples:

Text Only
hsw_bounds(35, 45, (1-0.27)*0.8, 75, 31, (1-0.27)*0.2, 2.2, 0) # porosity is 0.27
>>> (array([26.43276985]), array([7.07415429]), array([24.61588052]), array([0.]))
Source code in digirock/models/_hsw.py
Python
def hsw_bounds(
    bulk_mod: NDArrayOrFloat,
    shear_mod: NDArrayOrFloat,
    frac: NDArrayOrFloat,
    *argv: NDArrayOrFloat,
) -> Tuple[NDArrayOrFloat, NDArrayOrFloat, NDArrayOrFloat, NDArrayOrFloat]:
    """Hashin-Shtrikman-Walpole bounds

    $$
    K^{HS^+} = \\Lambda(\\mu_{max})
    $$

    Args:
        bulk_mod: first component bulk modulus
        shear_mod: first component shear modulus
        frac: first component volume fraction
        *argv: additional triplets of bulk_mod, shear_mod and frac for additional components

    Returns:
        K_hsp, K_hsm, Mu_hsp, Mu_hsm


    Refs:
        Rock Physics Handbook p171-172

    This example comes form the RPH p172

    Examples:
    ```
    hsw_bounds(35, 45, (1-0.27)*0.8, 75, 31, (1-0.27)*0.2, 2.2, 0) # porosity is 0.27
    >>> (array([26.43276985]), array([7.07415429]), array([24.61588052]), array([0.]))
    ```
    """
    args = _process_vfrac(*(bulk_mod, shear_mod, frac) + argv, i=2)
    shear_min = _hsw_get_minmax(
        *(v for i, v in enumerate(args) if i % 3 != 0), minmax="min"
    )
    shear_max = _hsw_get_minmax(
        *(v for i, v in enumerate(args) if i % 3 != 0), minmax="max"
    )
    bulk_min = _hsw_get_minmax(
        *(v for i, v in enumerate(args) if i % 3 != 1), minmax="min"
    )
    bulk_max = _hsw_get_minmax(
        *(v for i, v in enumerate(args) if i % 3 != 1), minmax="max"
    )
    k_hsp = _hsw_bulk_modulus_avg(
        shear_max, *(v for i, v in enumerate(args) if i % 3 != 1)
    )
    k_hsm = _hsw_bulk_modulus_avg(
        shear_min, *(v for i, v in enumerate(args) if i % 3 != 1)
    )
    zeta_hsp = _hsw_zeta(bulk_max, shear_max)
    zeta_hsm = _hsw_zeta(bulk_min, shear_min)
    mu_hsp = _hsw_shear_modulus_avg(
        zeta_hsp, *(v for i, v in enumerate(args) if i % 3 != 0)
    )
    mu_hsm = _hsw_shear_modulus_avg(
        zeta_hsm, *(v for i, v in enumerate(args) if i % 3 != 0)
    )
    to_shp = check_broadcastable(
        **{
            "k_hsp": k_hsp,
            "k_hsm": k_hsm,
            "mu_hsp": mu_hsm,
            "mu_hsm": mu_hsm,
        }
    )

    return tuple(np.broadcast_to(ar, to_shp) for ar in (k_hsp, k_hsm, mu_hsp, mu_hsm))

Cemented Sand Model

Cemented Sand Model

Refs

Mavko et. al. (2009) The Rock Physics Handbook, Cambridge Press pp 255-257

dryframe_cemented_sand(k_sand, mu_sand, k_cem, mu_cem, phi0, phi, ncontacts=9, alpha='scheme1')

Dryframe moduli for Cemented-Sand model.

Assumes sand grains are in contact and cement fills the porespace between grains.

Scheme 1: Cement preferentially fills spaces around contacts increasing shear strength. Scheme 2: Cement covers grains evenly but maintains contacts between sand grains. alpha can also be NDArrayOrFloat for a custom scheme.

Parameters:

Name Type Description Default
k_sand Union[numpy.ndarray, float]

Bulk modulus of sand grains

 required
mu_sand Union[numpy.ndarray, float]

Shear modulus of sand grains

 required
k_cem Union[numpy.ndarray, float]

Bulk modulus of cement

 required
mu_cem Union[numpy.ndarray, float]

Shear modulus of cement

 required
phi0 Union[numpy.ndarray, float]

Material porosity without any cement. For perfect spherical grains phi0 = 0.36

 required
phi Union[numpy.ndarray, float]

Material porosity after substituting in cement to porespace

 required
ncontacts int

Number of contacts between grains

 9
alpha Union[str, numpy.ndarray, float]

The cement fill scheme to use, one of ['scheme1', 'scheme2', NDArrayOrFloat]

 'scheme1'

Returns:

Type Description
Tuple[Union[numpy.ndarray, float], Union[numpy.ndarray, float]]

Bulk and shear modului of cemented-sand.
Source code in digirock/models/_cemented_sand.py
Python
def dryframe_cemented_sand(
    k_sand: NDArrayOrFloat,
    mu_sand: NDArrayOrFloat,
    k_cem: NDArrayOrFloat,
    mu_cem: NDArrayOrFloat,
    phi0: NDArrayOrFloat,
    phi: NDArrayOrFloat,
    ncontacts: int = 9,
    alpha: Union[str, NDArrayOrFloat] = "scheme1",
) -> Tuple[NDArrayOrFloat, NDArrayOrFloat]:
    """Dryframe moduli for Cemented-Sand model.

    Assumes sand grains are in contact and cement fills the porespace between grains.

    Scheme 1: Cement preferentially fills spaces around contacts increasing shear strength.
    Scheme 2: Cement covers grains evenly but maintains contacts between sand grains.
    alpha can also be `NDArrayOrFloat` for a custom scheme.

    Args:
        k_sand: Bulk modulus of sand grains
        mu_sand: Shear modulus of sand grains
        k_cem: Bulk modulus of cement
        mu_cem: Shear modulus of cement
        phi0: Material porosity without any cement. For perfect spherical grains phi0 = 0.36
        phi: Material porosity after substituting in cement to porespace
        ncontacts: Number of contacts between grains
        alpha: The cement fill scheme to use, one of ['scheme1', 'scheme2', `NDArrayOrFloat`]

    Returns:
        Bulk and shear modului of cemented-sand.
    """
    pois_sand = poisson_ratio(k_sand, mu_sand)
    pois_cem = poisson_ratio(k_cem, mu_cem)
    alpha = _cemented_sand_alpha(phi0, phi, ncontacts, alpha)

    lam_tau = mu_cem / (pi * mu_sand)
    lam_n = (2 * mu_cem * (1 - pois_sand) * (1 - pois_cem)) / (
        pi * mu_sand * (1 - 2 * pois_cem)
    )
    pois_sand2 = power(pois_sand, 2)
    c_tau = (
        10e-4
        * (9.654 * pois_sand2 + 4.945 * pois_sand + 3.1)
        * power(lam_tau, 0.01867 * pois_sand2 + 0.4011 * pois_sand - 1.8186)
    )
    b_tau = (0.0573 * pois_sand2 + 0.0937 * pois_sand + 0.202) * power(
        lam_tau, 0.0274 * pois_sand2 + 0.0529 * pois_sand - 0.8765
    )
    a_tau = (
        -0.01
        * (2.26 * pois_sand2 + 2.07 * pois_sand + 2.3)
        * power(lam_tau, 0.079 * pois_sand2 + 0.1754 * pois_sand - 1.342)
    )
    s_tau = a_tau * power(alpha, 2) + b_tau * alpha + c_tau
    c_n = 0.00024649 * power(lam_n, -1.9864)
    b_n = 0.20405 * power(lam_n, -0.89008)
    a_n = -0.024153 * power(lam_n, -1.3646)
    s_n = a_n * power(alpha, 2) + b_n * alpha + c_n

    keff = (ncontacts * k_cem * s_n / 6) * (1 - phi0)
    return (keff, 3 * keff / 5 + (3 * ncontacts * mu_cem * s_tau / 20.0) * (1 - phi0))