SOSAT package

Module contents

The State of Stress Analysis Tool (SOSAT) is a Python package that helps analyze the state of stress in the subsurface using various types of commonly available characterization data such as well logs, well test data such as leakoff and minifrac tests, regional geologic information, and constraints on the state of stress imposed by the existence of faults and fractures with limited frictional shear strength. It employs a Bayesian approach to integrate these data into a probability density function for the principal stress components.

Submodules

SOSAT.stress_state module

class SOSAT.stress_state.StressState(depth, avg_overburden_density, pore_pressure, depth_unit='m', density_unit='kg/m^3', pressure_unit='MPa', min_stress_ratio=0.4, max_stress_ratio=3.25, nbins=200, stress_unit='MPa')[source]

Bases: object

A class to contain all data necessary to define the probability distribution for all possible stress states at a given point in the subsurface.

depthfloat: The true vertical depth of the point being analyzed
vertical_stressfloat: The total vertical stress, currently taken as deterministic
pore_pressurefloat: The pore pressure, currently taken as deterministic
shmin_gridnumpy MaskedArray object: A 2D array containing the value of the minimum horizontal stress for each stress state considered. Values that are not admissible, such as where the minimum horizontal stress would be greater than the maximum horizontal stress are masked
shmax_gridnumpy MaskedArray object: A 2D array containing the value of the mmaximum horizontal stress for each stress state considered. Values that are not admissible, such as where the minimum horizontal stress would be greater than the maximum horizontal stress are masked
stress_unitstr: The unit used for stress
depth_unitstr: The unit used for vertical depth

depthfloat: the true vertical depth of the point being analyzed
average_overburden_densityfloat: average mass density of all overlying formations
pore_pressurefloat: formation pore pressure
depth_unitstr, optional: unit of measurement for depth, see list of units in pint package documentation
density_unitstr, optional: unit of measurement for mass density, see list of units in pintpackage documentation
pressure_unitstr, optional: unit of measurement for pressure, see list of units in pint package documentation
min_stress_ratiofloat, optional: Minimum stress included in the analysis expressed as a fraction of the vertical stress. This is not intended as a way to express a prior probability but is instead meant as a convenient way to set the bounds of stress considered in the analysis. It does effectively truncate the stress distribution at this value and so should be chose well outside of the zone of significant probability density. Default value is 0.4
max_stress_ratiofloat: Maximum stress included in the analysis expressed as a fraction of the vertical stress. This is not intended as a way to express a prior probability but is instead meant as a convenient way to set the bounds of stress considered in the analysis. It does effectively truncate the stress distribution at this value and so should be chose well outside of the zone of significant probability density. Default value is 3.25
nbinsint: number of bins to use for the horizontal principal stresses. The same number is used for both horizontal principal stresses. Default value is 200.

This class applies the Bayesian approach to quantifying uncertainty in the state of stress at a single point in the subsurface that is outlined in [1]_. By “single point” here we mean any volume of the subsurface that the user wishes to treat as having a homogeneous stress state the. Depending on the application this could be a volume with length scale of a log scale (tens of centimeters) or in a regional-scale analysis could be a volume with length scale of kilometers.

Currently the vertical stress and pore pressure are taken as deterministically known, so that uncertainties in these values are not reflected in the uncertainty in the the horizontal principal stresses. This may change in a future release, but in many applications is not a big limitation because there is generally a smaller uncertainty in the vertical stress and pore pressure compared to the horizontal principal stresses.

Currently this class does not contain stress orientation information, though this too is planned to change in a future release. Though this would only add significant value when methods to add different StressState objects are needed and implemented. This will be required, for example, to interpolate between different points where a StressState object has been parameterized.

[1] Burghardt, J. “Geomechanical Risk Assessment for Subsurface Fluid Disposal Operations,” Rock Mech Rock Eng 51, 2265–2288 (2018) https://doi.org/10.1007/s00603-018-1409-1

To compute and plot the posterior distribution at a point with a frictional faulting constraint, you would do the following:

>>> from SOSAT import StressState
>>> from SOSAT.constraints import FaultConstraint
>>> ss = StressState(1.0,
                     2.5,
                     0.3,
                     depth_unit='km',
                     density_unit='g/cm^3',
                     pressure_unit='MPa')

>>> sigv = ss.vertical_stress
>>> fc = FaultConstraint()
>>> ss.add_constraint(fc)
>>> fig = ss.plot_posterior()
>>> plt.savefig("fault_constraint_posterior.png")

A_phi_calculate()[source]

Computes the scalar A_phi value for each stress state included in the class. The scalar regime parameter varies from 0 to 1. The scalar measure of the faulting regime is defined in Simpson (1997):

phi = (sigma2-sigma3)/ (sigma1-sigma3) n = 0, 1, and 2 respectively for normal, strike-slip, and reverse faulting regime respectively. A_phi = (n+0.5) + (-1)^n*(phi-0.5)

Using this definition, A_phi value ranges from 0 to 3, where [0,1] corresponds to normal faulting regime, [1,2] corresponds to strike-slip faulting regime, and [2,3] corresponds to reverse faulting regime.

2D Numpy MaskedArray with dtype=float: The array is square with dimensions nbins

add_constraint(constraint)[source]

Method to add a constraint to the stress state probability distribution. Constraints will be applied in the order that they are added

constraintConstraint object: The object containing the information for the constraint, must be a class in the SOSAT.constraints submodule

evaluate_posterior()[source]

Method to evaluate the posterior joint probability density given the constraints that have been added. If no constraints have been added the prior distribution will be returned.

2D Numpy MaskedArray with dtype=float: The array is square with dimensions nbins, where each entry in the array contains the probability density for the stress state defined by shmin_grid and shmax_grid

The probability density will be properly normalized such that the the values of all bins in the distribution will sum to one. The prior distribution is homogeneous for all compressive stress states considered in the distribution. The stress states considered in the distribution are governed by the min_stress_ratio and max_stress_ratio passed into the constructor of this class.

The posterior distribution is computed by recursively applying Bayes’ law using the likelihood function provided by each constraint object that has been added to this class.

get_shmax_confidence_intervals(confidence)[source]: return an upper and lower bound within a specified confidence interval. Confidence interval should be specified as a fraction so that, for example, a 95% confidence interval is expressed as confidence=0.95

get_shmax_marginal()[source]: get the marginal probability distribution for the maximum principal stress

get_shmax_marginal_cdf()[source]: get the marginal cumulative probability function for the maximum horizontal stress

get_shmin_confidence_intervals(confidence)[source]: return an upper and lower bound within a specified confidence interval. Confidence interval should be specified as a fraction so that, for example, a 95% confidence interval is expressed as confidence=0.95

get_shmin_marginal()[source]: get the marginal probability distribution for the minimum principal stress

get_shmin_marginal_cdf()[source]: get the marginal cumulative probability function for the minimum horizontal stress

plot_posterior(figwidth=5.0, figheight=3.5, contour_levels=5, cmap=<matplotlib.colors.LinearSegmentedColormap object>)[source]

Makes a contour plot of the joint probability distribution of the maximum and minimum horizontal stresses

figwidthfloat, optional: the width of the figure in inches, defaults to 5 inches
figheightfloat, optional: the height of the figure in inches, defaults to 5 inches
contour_levelsint, optional: the number of contour levels desired in the plot
cmapcolormap object, optional: a matplotlib colormap object to use to display the probability density. Default is matplotlib.pyplot.cm.Greys

matplotlib.pyplot.Figure object

The plot is generated with matplotlib.pyplot.contourf

regime()[source]

Computes the scalar regime parameters for each stress state included in the class. The scalar regime parameter varies from negative one to one. It is defined the vector space where each each stress state is represented by a vector whose head lies at the coordinate (shmax,shmin) and whose tail lies at the spherical stress state where

shmax = shmin = vertical stress

The scalar measure of the faulting regime is defined by the dot product of the vector for the given stress state and the unit vector give by

shmax = shmin = - 1/sqrt(2)

Using this definition values of the regime parameter between -1 and -sqrt(2)/2 correspond to thrust faulting states, values between -sqrt(2)/2 and sqrt(2)/2 correspond to strike- slip states, and values between sqrt(2)/2 and +1 correspond to normal faulting states.

2D Numpy MaskedArray with dtype=float: The array is square with dimensions nbins

SOSAT.stress_state.fmt(x, pos)[source]: A utility function to improve the formatting of plot labels

SOSAT.stress_state.units = <pint.registry.UnitRegistry object>: The stress_state module contains classes and functions that are used to represent and plot the joint probability distribution for the minimum and maximum horizontal stresses.

SOSAT.constraints

SOSAT.constraints.fault_constraint

class SOSAT.constraints.fault_constraint.FaultConstraint(friction_dist=None)[source]

Bases: SOSAT.constraints.constraint_base.StressConstraint

A class used to constrain the state of stress by the existence of frictional fault and fractures.

friction_distsubclass of scipy.stats.rv_continuous: The probability density function for the fault friction coefficient

friction_distsubclass of scipy.stats.rv_continuous, optional: a probability distribution function from the fault friction coefficient. If not provided the default is a lognormal distribution with s=0.15 and scale=0.7 passed into scipy.stats.lognorm

loglikelihood(ss)[source]

Computes the likelihood of each stress state

ss: SOSAT.StressState object: StressState object containing the stress states over which the likelihood is to be evaluated

Numpy MaskedArray: The returned object is a Numpy MaskedArray containing the likelihood for each stress ss. The returned array is masked identically to ss.shmin_grid

SOSAT.constraints.fault_constraint.critical_friction(sig1, sig3, pp)[source]

A function that computes the critical friction coefficient that would induce slip for the given combination of minimum and maximum principal stress and pore pressure

sig1float or array_like: The maximum (most compressive) principal stress. If sig1 is array_like then sig3 must be a scalar of type float
sig3float or array_like: The minimum (least compressive) principal stress. If sig3 is array_like then sig1 must be a scalar of type float
ppfloat: The pore pressure

mu_cfloat or array_like: The critical friction coefficient. If either sig1 or sig3 are array_like then the result will be an array containing the critical friction coefficient for each entry in the array.

SOSAT.constraints.faulting_regime_constraint

class SOSAT.constraints.faulting_regime_constraint.FaultingRegimeConstraint(regime_dist)[source]

Bases: SOSAT.constraints.constraint_base.StressConstraint

A class that constrains the state of stress based on knowledge or assumptions about the relative probability of the three Andersonian faulting regimes.

regime_distsubclass of scipy.stats.rv_continuous: The probability density function for the scalar regime parameter. See StressState.regime for details on the definition of the scalar regime parameter. Usually the and instance of SU is provided to apply a uniform probability to each Andersonian faulting regime with logistic functions to smooth the transition between each regime. However, any subclass of scipy.stats.rv_continuous defined from -1.0 to 1.0 can be used.

loglikelihood(ss)[source]

Abstract method that requires all derived classes to implement a funcction that computes the loglikelihood of the constraint for a the stress states passed in through the StressState parameter

ss: SOSAT.StressState object: StressState object containing the stress states over which the likelihood is to be evaluated

Numpy MaskedArray: The returned object is a Numpy MaskedArray containing the log likelihood for each stress ss. The returned array is masked identically to ss.shmin_grid

class SOSAT.constraints.faulting_regime_constraint.SU_gen(momtype=1, a=None, b=None, xtol=1e-14, badvalue=None, name=None, longname=None, shapes=None, extradoc=None, seed=None)[source]

Bases: scipy.stats._distn_infrastructure.rv_continuous

A class used to generate smoothed semi-uniform probability distributions for the scalar regime parameter. Most commonly this is used to define a probability distribution for the three Andersonian faulting regimes. A logistic function is used to smooth the transition between the different faulting regimes.

w_NFfloat: the weight assigned to stress states for which the scalar regime parameter is less than k2; with the default value of k2 this corresponds to normal faulting stress states; the absolute value only matters relative to the weight assigned to the other two faulting regimes
w_SSfloat: the weight assigned to stress states for which the scalar regime parameter is between k1 and k2; with the default values of k1 and k2 this corresponds to strike-slip faulting stress states; the absolute value only matters relative to the weight assigned to the other two faulting regimes
w_TFfloat: the weight assigned to stress states for which the scalar regime parameter is greater than k1; with the default value of k1 this corresponds to thrust faulting stress states; the absolute value only matters relative to the weight assigned to the other two faulting regimes
theta1float, optional: the first cutoff angle; defaults to a value of sqrt(2)/2, which is the value of the scalar regime parameter that marks the transition from thrust faulting to strike-slip faulting regimes
k1float, optional: a parameter controlling the width of the sigmoid transition at theta1. Larger numbers yield a more gradual transition. Default value is 300
theta2float: the second cutoff angle; default value is -sqrt(2)/2, which is the value of the scalar regime parameter that makrs the transition from strike-slip faulting to normal faulting.
k2foat, opptional: a parameter controlling the width of the sigmoid transition at theta2. Larger numbers yield a more gradual transition. Default value is 300

The methods of scipy.stats.rv_continuous that are overridden are private methods. The corresponding public methods check their arguments before calling the private methods.

SOSAT.constraints.faulting_regime_constraint.int_logist(k, xo)[source]

SOSAT.constraints.faulting_regime_constraint.int_logist_a(a, k, xo)[source]

SOSAT.constraints.faulting_regime_constraint.logistic(x, k, xo)[source]

SOSAT.constraints.breakout_constraint

class SOSAT.constraints.breakout_constraint.BreakoutConstraint(breakout_exists, UCS_dist, rock_friction_angle_dist, rock_friction_angle_units, mud_pressure_dist, mud_temperature_dist, formation_temperature, YM, PR, CTE, pressure_unit='Pa')[source]

Bases: SOSAT.constraints.constraint_base.StressConstraint

A class used to constrain the stress state by the existence or non existence of borehole breakouts at the location being analyzed. If breakouts exists that generally indicates that the maximum horizontal stress is much larger than the minimum horizontal stress, if not breakout exists that puts a limit on how large the maximum horizontal stress could be; the magnitude of that limit depends on the strength of the rock, drilling mud weight, temperature, etc. Significant uncertainty in these parameters, especially the rock strength and mud weight makes the constraint provided by a breakout analysis much weaker. This class seeks to account for these uncertainties in the posterior stress distribution

No public attributes

breakout_existsbool: Indication whether or not breakouts exist
UCS_distsubclass of scipy.stats.rv_continuous: A probability distribution for the unconfined compressive strength. Shoudl use the same pressure unit as is used by the mud pressure and Young’s modulus, but this can be any unit as specified though the optional pressure_unit parameter; this shoudl be an estimate for the minimum UCS in the zone of interest since that is what would govern the formation of breakouts, not the average or representative value
rock_friction_angle_distsubclass of scipy.stats.rv_continuous: A probabiilty distribution for the rock friction angle.
rock_friction_angle_unitsstr: The unit used for the rock friction angle. Should be an angular unit recognized by pint.UnitRegistry (i.e. ‘deg’ for degrees, ‘radians’ for radians
mud_pressure_distsubclass of scipy.stats.rv_continuous: The probability distribution for the minimum mud pressure experienced by the relevant section of borehole from the time that the well was drilled until the log used to identify the presence or absence of breakouts was run; mud pressure should be specified in the same pressure unit as is used for UCS and Young’s modulus, but this can be any unit as specified though the optional pressure_unit parameter, which defaults to ‘Pa’; conversion from mud weight must be performed by the user of this class
mud_temperature_distsubclass of scipy.stats.rv_continuous: The probability distribution for the maximum mud temperature; the maximum value is of interest rather than the average value since the formation of a breakout is governed by the maximum value only
formation_temperaturefloat: Formation temperature, which is taken as deterministic since it is usually not highly uncertain
YMfloat: Formation Young’s Modulus, which is taken as deterministic since the formation of breakouts is only weakly dependent on this parameter; should be specified in the same pressure unit as is used for mud pressure and Young’s modulus, but this can be any unit as specified though the optional pressure_unit parameter, which defaults to ‘Pa’
PRfloat: Formation Poisson’s Ratio, which is taken as deterministic since the formation of breakouts is only weakly dependent on this parameter
CTEfloat: Formation coefficient of thermal expansion, which is taken as deterministic since the formation of breakouts is only weakly dependent on this parameter
pressure_unitstr, optional: The unit used for UCS and Young’s modulus; should be a unit recognized by pint.UnitRegistry; defaults to ‘Pa’

This class uses a Mohr Coulomb (MC) failure criterion to evaluate the probability of a breakout occuring. The Mohr Coulomb criterion uses two parameters, most often a friction angle and cohesion. In this class the unconfined compressive strength (UCS) and friction angle were chosen instead of the cohesion and friction angle. This is because the UCS is a more commonly known and intuitive property, and the cohesion MC parameter can be calculated given the UCS and friction angle.

While this class allows users to use any probability distribution that derives from the scipy.stats.rv_continuous class, users are cautioned against using any distribution that has finite probability density for negative parameter values, since negative strength values are not physically meaningful. Therefore, lognormal distributions are more appropriate than a normal distribution, for example.

loglikelihood(ss)[source]

Computes the loglikelihood of each stress state given the presence or absence of breakouts, formation properties, and mud properties specified

ss: SOSAT.StressState object: StressState object containing the stress states over which the likelihood is to be evaluated

Numpy MaskedArray: The returned object is a Numpy MaskedArray containing the likelihood for each stress ss. The returned array is masked identically to ss.shmin_grid

SOSAT.constraints.breakout_constraint.MCfail(sig1, sig3, phi, Co)[source]

Mohr Coulomb failure criterion

sig1float, array_like: maximum principal stress
sig3float, array_like: minimum principal stress
phifloat: friction angle
Cofloat: cohesion

ffloat or array_like: The value of the failure function, which is defined so that negative values indicate no failure and positive values indicate failure, a zero value is at incipient failure

SOSAT.constraints.DITF_constraint

class SOSAT.constraints.DITF_constraint.DITFConstraint(DITF_exists, mud_pressure_dist, mud_temperature_dist, tensile_strength_dist, formation_temperature, YM, PR, CTE, pressure_unit='Pa')[source]

Bases: SOSAT.constraints.constraint_base.StressConstraint

A class used to constrain the stress state by the existence or non existence of drilling-induced tensile fractures (DITF) at the location being analyzed. Depending on the mud and formation temperatures, mud weights, and rock strength and whether or not significant mud losses were observed, if DITF’s exist it generally indicates that the maximum horizontal stress is much larger than the minimum principal stress.

No public attributes

DITF_existsbool: Indication whether or not DITF exist
mud_pressure_distsubclass of scipy.stats.rv_continuous: The probability distribution for the maximum mud pressure experienced by the relevant section of borehole from the time that the well was drilled until the log used to identify the presence or absence of breakouts was run; mud pressure should be specified in the same pressure unit as is used for UCS and Young’s modulus, but this can be any unit as specified though the optional pressure_unit parameter, which defaults to ‘Pa’; conversion from mud weight must be performed by the user of this class
mud_temperature_distsubclass of scipy.stats.rv_continuous: The probability distribution for the minimum mud temperature; the minimum value is of interest rather than the average value since the formation of a DITF is governed by the minimum value only
tensile_strength_distsubclass of scipy.stats.rv_continuous: The probability distribution for minimum the minimum tensile strength in the zone being analyzed. DITFs will form at the weakest portion of the well for a given stress state, so whether they form or not is dependent on the minimum tensile strength rather than an average representative value
formation_temperaturefloat: Formation temperature, which is taken as deterministic since it is usually not highly uncertain
YMfloat: Formation Young’s Modulus, which is taken as deterministic since the formation of DITF is only weakly dependent on this parameter; should be specified in the same pressure unit as is used for mud pressure and Young’s modulus, but this can be any unit as specified though the optional pressure_unit parameter, which defaults to ‘Pa’
PRfloat: Formation Poisson’s Ratio, which is taken as deterministic since the formation of DITFs is only weakly dependent on this parameter
CTEfloat: Formation coefficient of thermal expansion, which is taken as deterministic since the formation of DITF is only weakly dependent on this parameter
pressure_unitstr, optional: The unit used for UCS and Young’s modulus; should be a unit recognized by pint.UnitRegistry; defaults to ‘Pa’

While this class allows users to use any probability distribution that derives from the scipy.stats.rv_continuous class for the mud temperature, pressure, and formation tensile strength, users are cautioned against using any distribution that has finite probability density for negative parameter values, since negative strength values are not physically meaningful. Therefore, lognormal distributions are more appropriate than a normal distribution, for example.

loglikelihood(ss)[source]

Computes the likelihood of each stress state given the presence or absence of DITFs, formation and mud properties specified.

ss: SOSAT.StressState object: StressState object containing the stress states over which the likelihood is to be evaluated

Numpy MaskedArray: The returned object is a Numpy MaskedArray containing the likelihood for each stress ss. The returned array is masked identically to ss.shmin_grid

SOSAT.constraints.stress_measurement_constraint

class SOSAT.constraints.stress_measurement_constraint.StressMeasurement(shmin_dist)[source]

Bases: SOSAT.constraints.constraint_base.StressConstraint

A class to incorporate mini-frac or extended leakoff test (XLOT) data to constrain the state of stress

No public attributes

shmin_distsubclass of scipy.stats.rv_continuous: The probability distribution for the measured minimum horizontal stress

This assumes that the minimum principal stress, which is measured by mini-frac or XLOT tests is horizontal. This makes this class only applicable to normal faulting or strike-slip stress environments

loglikelihood(ss)[source]

Compute the likelihood of each stress state in ss based on the measured stress

ssSOSAT.StressState object: The stress states on which to evaluate the likelihood

SOSAT.samplers.rejection_sampler

class SOSAT.samplers.rejection_sampler.RejectionSampler(SS)[source]

Bases: object

A class used to generate samples of the posterior stress distribution using a rejection sampling approach

SS : SOSAT.StressState object

GenerateSamples(Nsamples)[source]

Method to evaluate samples

Nsamples: int: Number of samples to generate

shmin, shmax, sv: arrays containing samples of the three: principal stress

SOSAT.samplers.rejection_sampler.sigHcritTF(sigv, pp, mu)[source]

SOSAT.samplers.rejection_sampler.sighcritNF(sigv, pp, mu)[source]

SOSAT.cli module

Console script for SOSAT. Not yet implemented