Anomaly quantification

reheatfunq.anomaly

The reheatfunq.anomaly module contains functionality to analyze the strength of heat flow anomalies using the GammaConjugatePrior model of regional aggregate heat flow distributions. The module contains the workhorse HeatFlowAnomalyPosterior for Bayesian heat flow anomaly strength quantification and AnomalyLS1980 class to model a fault-generated conductive heat flow anomaly [LS1980]. The class AnomalyNearestNeighbor can be used to include the results of an external heat flow anomaly modeling (finite elements etc.) in the REHEATFUNQ analysis. The workflow for anomaly quantification using REHEATFUNQ consists of the following steps:

1. Define the \(d_\mathrm{min}\) (e.g. \({d_\mathrm{min} = 20\,\mathrm{km}}\))

2. Define the conjugate prior to use. Obtain a GammaConjugatePrior instance (e.g. using the REHEATFUNQ default from default_prior()).

3. Model the fault-generated heat flow anomaly. So far, the AnomalyLS1980 and AnomalyNearestNeighbor are available for this purpose.

4. Compute the marginal posterior in \(P_H\) using the HeatFlowAnomalyPosterior class, which takes into consideration the bootstrapped updating of the gamma conjugate prior over the set of \(d_\mathrm{min}\)-conforming subsets of the heat flow data.

Vertical Strike-Slip Fault

Exemplarily, the following code summarizes the analysis using a heat flow anomaly for a vertical strike-slip fault [LS1980]. First, we generate some toy heat flow data following a gamma distribution. We use the same heat flow data as in the reheatfunq.regional example:

import numpy as np
from reheatfunq.anomaly import AnomalyLS1980
rng = np.random.default_rng(123920)
alpha = 53.3
qu = rng.gamma(alpha, size=15)
x = 100e3 * (rng.random(15) - 0.5)
y = 100e3 * (rng.random(15) - 0.5)

Generate an obliquely striking vertical strike slip fault and the corresponding conductive heat flow anomaly for a linearly increasing heat production with depth [LS1980]:

fault_trace = np.array([(-20e3, -100e3), (20e3, 100e3)])
anomaly = AnomalyLS1980(fault_trace, 14e3)
xy = np.stack((x,y), axis=1)

The data and anomaly look like this (dashed black lines indicate the contours of the heat flow anomaly \(c_i=\Delta q_i / P_H\) and the blue line shows the fault trace):

Now compute three sets of heat flow data superposed by heat flow anomalies of \(90\,\mathrm{MW}\), \(150\,\mathrm{MW}\) and \(300\,\mathrm{MW}\) power:

dq = anomaly(xy)
q1 = qu +  90e6 * dq * 1e3
q2 = qu + 150e6 * dq * 1e3
q3 = qu + 300e6 * dq * 1e3

Now compute the marginalized posterior distribution of the heat-generating power \(P_H\) for the data superposed with the three anomalies:

from reheatfunq.anomaly import HeatFlowAnomalyPosterior
from reheatfunq import default_prior
gcp = default_prior()
post1 = HeatFlowAnomalyPosterior(q1, x, y, anomaly, gcp)
post2 = HeatFlowAnomalyPosterior(q2, x, y, anomaly, gcp)
post3 = HeatFlowAnomalyPosterior(q3, x, y, anomaly, gcp)

P_H = np.linspace(0, 600e6, 200)
pdf1 = post1.pdf(P_H)
pdf2 = post2.pdf(P_H)
pdf3 = post3.pdf(P_H)

The vertical dashed lines indicate the true anomaly powers.

A detailed use of the anomaly quantification can be found in the Jupyter notebook jupyter/REHEATFUNQ/06-Heat-Flow-Analysis.ipynb.

Custom Heat Flow Anomaly

The next example shows how to use REHEATFUNQ with a heat flow anomaly that has been computed using external code. This is done using the AnomalyNearestNeighbor class, and the heat flow anomaly should be specified in terms of the factors \(c_i\) (i.e. \(\Delta q_i / P_H\)). For this purpose, we generate a Gauss-shaped heat flow anomaly that leads to an additional heat flow of \(68.3\,\mathrm{mW\,m}^{-2}\) at its center when fed by \(10\,\mathrm{MW}\) of heat-generating power:

def anomaly_ci(x,y):
    return 68.3e-3 / 10e6 * np.exp(-(x**2 + y**2))

Note here that the \(c_i\) should be specified in basic SI units.

We generate some new data superposed by a \(10\,\mathrm{MW}\) heat flow anomaly:

N = 20
rng = np.random.default_rng(123329773)
xy = 3 * rng.random((N,2)) - 1.5
q0 = 0.39 * rng.gamma(175.2, size=N)
c_i = anomaly_ci(*xy.T)
q = q0 + 1e3 * 10e6 * c_i

To adjust the analysis to a custom set of heat flow anomaly factors \(c_i\) in \(\mathrm{m}^{-2}\), it is sufficient to replace the line c_i = anomaly_ci(*xy.T) with whatever code computes or loads the factors. The shape of c_i should be compatible to (N,).

The point set and anomaly generated by the above code should look like this:

With the anomaly factors \(c_i\) evaluated at the data locations, we can perform the REHEATFUNQ anomaly analysis using the AnomalyNearestNeighbor class:

from reheatfunq.anomaly import AnomalyNearestNeighbor
xyc = np.stack((*xy.T, c_i), axis=1)
ann = AnomalyNearestNeighbor(xyc)
hfap = HeatFlowAnomalyPosterior(q, *xy.T, ann, gcp, dmin=0.0)

The analysis recovers the anomaly strength:

class HeatFlowAnomalyPosterior(q, x, y, anomaly, gcp, dmin=20000.0, n_bootstrap='auto', heat_flow_unit='mW/m²', rng=127, rtol=1e-08, pdf_algorithm='barycentric_lagrange', bli_max_refinements=7, precision='long double')

This class evaluates the posterior probability of the strength of a heat flow anomaly, expressed by the frictional power \(P_H\) on the fault, using the REHEATFUNQ model of regional heat flow and a set of regional heat flow data.

Parameters:

• q (array_like) – The heat flow data of shape (N,).

• x (array_like) – The \(x\) locations of the heat flow data. Also shape (N,).

• y (array_like) – The \(y\) locations of the heat flow data. Also shape (N,).

• anomaly (reheatfunq.anomaly.anomaly.Anomaly | list[reheatfunq.anomaly.anomaly.Anomaly] | list[tuple[float,reheatfunq.anomaly.anomaly.Anomaly]]) –

  The model of the heat flow anomaly that can be evaluated at the data locations. There are three ways to specify this parameter:

  1. A single Anomaly instance.

  2. A list of Anomaly instances.

  3. A list of tuples (float,Anomaly).

  In the second and third case, the Anomaly instances are interpreted as a discretization of the heat transport uncertainty, that is, they should span the space of possible heat transport solution given the knowledge of the problem at hand. The third case lists float as prior propabilities for each of the anomaly. This way, one can specify a probability distribution of the heat transport solutions.

• gcp (reheatfunq.regional.GammaConjugatePrior | tuple) – The prior for the regional aggregate heat flow distribution.

• dmin (float, optional) – The minimum distance between data points (in m). If data points closer than this distance exist in the heat flow data, they are not considered independent and are alternated in the bootstrap.

• n_bootstrap (int, optional) – The number of permuted heat flow data sets to generate. If no pair of data points is closer than the minimum distance \(d_\mathrm{min}\), this parameter has no effect.

• heat_flow_unit ('mW/m²' | 'W/m²', optional) – The unit in which the heat flow data q are given.

• rng (int | numpy.random.Generator, optional) – The random number generator to use or a reproducible seed.

• rtol (float, optional) – The relative tolerance to aim for in various places within the underlying numerics (quadrature and interpolation).

• pdf_algorithm ("explicit" | "barycentric_lagrange" | "adaptive_simpson", optional) –

  The algorithm to use for evaluating the PDF.

  1. "explicit": Explicitly calculate all requested points.

  2. "barycentric_lagrange": First create a barycentric Lagrange interpolator of the PDF which is thereafter used to evaluate the PDF

  3. "adaptive_simpson": Create the adaptive Simpson’s rule integrator using the explicitly evaluated PDF. Then use the integrator’s polynomials to evaluate the PDF thereafter.

  The chosen PDF evaluation algorithm will also be used when evaluating the cumulative distribution functions.

• bli_max_refinements (int, optional) – The maximum number of refinements of the barycentric Lagrange interpolator. The refinement will stop before if the precision goal is reached. The number of Chebyshev sampling points of the explicitly evaluated PDF grows base-2 exponentially with bli_max_refinements.

• precision ('double' | 'long double', optional) – The precision of the internal numerical computations. The higher the precision, the more likely it is to obtain a precise result for large data sets. The trade off is a longer run time. If the respective flags have been set at compile time, additional options 'float128' (GCC 128bit floating point), 'dec50' (boost 50-digit multiprecision), and 'dec100' (boost 100-digit multiprecision) are available.

cdf(P_H)

Evaluate the marginal posterior cumulative distribution of heat-generating power \(P_H\).

Parameters:

P_H (array_like) – The powers (in W) at which to evaluate the marginal posterior cumulative distribution.

Returns:

cdf – The marginal posterior cumulative distribution of heat-generating power \(P_H\) evaluated at the given P_H.

Return type:

numpy.typing.NDArray[numpy.float64]

pdf(P_H)

Evaluate the marginal posterior distribution in heat-generating power \(P_H\).

Parameters:

P_H (array_like) – The powers (in W) at which to evaluate the marginal posterior density.

Returns:

pdf – The marginal posterior probability density of heat-generating power \(P_H\) evaluated at the given P_H.

Return type:

numpy.typing.NDArray[numpy.float64]

tail(P_H)

Evaluate the posterior tail distribution (complementary cumulative distribution) of heat-generating power \(P_H\).

Parameters:

P_H (array_like) – The powers (in W) at which to evaluate the marginal posterior tail distribution.

Returns:

tail – The marginal posterior tail distribution of heat-generating power \(P_H\) evaluated at the given P_H.

Return type:

numpy.typing.NDArray[numpy.float64]

tail_quantiles(quantiles)

Compute posterior tail quantiles, that is, heat-generating powers \(P_H\) at which the complementary cumulative distribution of \(P_H\) has fallen to level \(x\).

Parameters:

quantiles (array_like) – The tail quantiles to compute.

Returns:

P_H – The heat-generating power \(P_H\) at which the posterior tail distribution evaluates to x.

Return type:

numpy.typing.NDArray[numpy.float64]

class Anomaly

Base class for all heat flow anomalies. This base class is useless on its own, it only provides the call signature for the underlying C++ implementation of the anomaly evaluation. The inheriting classes provide the C++ implementation and inherit the __call__ functionality.

Inheriting classes:

• AnomalyLS1980.

__call__(xy, P_H=1.0)

Evaluate the heat flow anomaly at a set of points for a given heat-generating power P_H.

Parameters:

• xy (double[:,:]) – Locations at which to evaluate the heat flow anomaly.

• P_H (float) – The heat-generating power (in W).

Returns:

The anomalous heat flow evaluated at the locations, \(\{\Delta q_i\} = \{c_i P_H \}\).

Return type:

numpy.ndarray

class AnomalyLS1980(const double[:, ::1] xy, double d)

Bases: Anomaly

A conductive heat flow anomaly generated by a vertical strike-slip fault whose heat production increases linearly with depth.

This model uses equation (A23b) of Lachenbruch & Sass [LS1980] which is an analytical solution for a straight, infinitely long vertical strike-slip fault in a homogeneous half space. For each queried point, the closest point on the actual, segmented fault is computed using CGAL. The distance between this fault trace point and the query point is plugged into equation (A23b) (see the REHEATFUNQ paper for further details).

The quality of this model’s approximation depends on the curvature of the fault and the homogeneity of heat conduction in the crustal volume of interest.

Parameters:

• xy (array_like) – Array of consecutive fault trace coordinates of shape (N,2). The second dimension iterates the coordinate tuple (x[i], y[i]).

• d (float) – Depth of the fault (in m).

__call__(xy, P_H=1.0)

See Anomaly.

length(self) → double

class AnomalyNearestNeighbor(const double[:, ::1] xyc)

Bases: Anomaly

A multi-purpose heat flow anomaly class with user-provided anomaly coefficients \(c_i\) using nearest neighbor interpolation to obtain the coefficients at the data locations.

This class can be used to provide arbitrary heat flow solutions (e.g. from numerical methods) to REHEATFUNQ. To do so, the coefficients should be provided at the locations of the data later analyzed.

Parameters:

xyc (array_like) – Array of point solutions \((x_i,y_i,c_i)\).

__call__(xy, P_H=1.0)

See Anomaly.

The use of this class is demonstrated in the quickstart Jupyter notebook jupyter/Custom-Anomaly.ipynb.