Climate change signal in well temperature logs (fatiando.geothermal.climsig)¶

Modeling and inversion of temperature residuals measured in wells due to temperature perturbations in the surface.

Perturbations can be of two kinds: abrupt or linear.

Forward modeling of these types of changes is done with functions:

• abrupt
• linear

Assumeing that the temperature perturbation was abrupt. The residual temperature at a depth \(z_i\) in the well at a time \(t\) after the perturbation is given by

where \(A\) is the amplitude of the perturbation, \(\lambda\) is the thermal diffusivity of the medium, and \(\mathrm{erf}\) is the error function.

For the case of a linear change, the temperature is

Given the temperature measured at different depths, we can invert for the amplitude and age of the change. The available inversion solvers are:

• SingleChange: inverts for the parameters of a single temperature change. Can use both abrupt and linear models.

class fatiando.geothermal.climsig.SingleChange(temp, zp, mode, diffus=31.5576)[source]¶

Bases: fatiando.inversion.misfit.Misfit

Invert the well temperature data for a single change in temperature.

The parameters of the change are its amplitude and age.

See the docstring of fatiando.geothermal.climsig for more information and examples.

Parameters:

• temp

  : array

  The temperature profile

• zp

  : array

  Depths along the profile

• mode

  : string

  The type of change: 'abrupt' for an abrupt change, 'linear' for a linear change.

• diffus

  : float

  Thermal diffusivity of the medium (in m^2/year)

Note

The recommended solver for this inverse problem is the Levemberg-Marquardt method. Since this is a non-linear problem, set the desired method and initial solution using the config method. See the example bellow.

Example with synthetic data:

>>> import numpy as np
>>> zp = np.arange(0, 100, 1)
>>> # For an ABRUPT change
>>> amp = 2
>>> age = 100 # Uses years to avoid overflows
>>> temp = abrupt(amp, age, zp)
>>> # Run the inversion for the amplitude and time
>>> # This is a non-linear problem, so use the Levemberg-Marquardt
>>> # algorithm with an initial estimate
>>> solver = SingleChange(temp, zp, mode='abrupt')
>>> _ = solver.config('levmarq', initial=[1, 1]).fit()
>>> amp_, age_ = solver.estimate_
>>> print("amp: {:.2f}  age: {:.2f}".format(amp_, age_))
amp: 2.00  age: 100.00
>>> # For a LINEAR change
>>> amp = 3.45
>>> age = 52.5
>>> temp = linear(amp, age, zp)
>>> solver = SingleChange(temp, zp, mode='linear')
>>> _ = solver.config('levmarq', initial=[1, 1]).fit()
>>> amp_, age_ = solver.estimate_
>>> print("amp: {:.2f}  age: {:.2f}".format(amp_, age_))
amp: 3.45  age: 52.50

Notes:

For abrupt changes, derivatives with respect to the amplitude and age are calculated using the formula

\[\frac{\partial T_i}{\partial A} = 1 - \mathrm{erf}\left( \frac{z_i}{\sqrt{4\lambda t}}\right)\]

and

\[\frac{\partial T_i}{\partial t} = \frac{A}{t\sqrt{\pi}} \left(\frac{z_i}{\sqrt{4\lambda t}}\right) \exp\left[-\left(\frac{z_i}{\sqrt{4\lambda t}}\right)^2\right]\]

respectively.

For linear changes, derivatives with respect to the age are calculated using a 2-point finite difference approximation. Derivatives with respect to amplitude are calculate using the formula

\[\frac{\partial T_i}{\partial A} = \left(1 + 2\frac{z_i^2}{4\lambda t}\right) \mathrm{erfc}\left(\frac{z_i}{\sqrt{4\lambda t}}\right) - \frac{2}{\sqrt{\pi}}\left(\frac{z_i}{\sqrt{4\lambda t}}\right) \mathrm{exp}\left(-\frac{z_i^2}{4\lambda t}\right)\]

config(method, **kwargs)¶

Configure the optimization method and its parameters.

This sets the method used by fit and the keyword arguments that are passed to it.

Parameters:

• method

  : string

  The optimization method. One of: 'linear', 'newton', 'levmarq', 'steepest', 'acor'

Other keyword arguments that can be passed are the ones allowed by each method.

Some methods have required arguments:

• newton, levmarq and steepest require the initial argument (an initial estimate for the gradient descent)
• acor requires the bounds argument (min/max values for the search space)

See the corresponding docstrings for more information:

• linear
• newton
• levmarq
• steepest
• acor

copy(deep=False)¶

Make a copy of me together with all the cached methods.

estimate_¶

A nicely formatted version of the estimate.

If the class implements a fmt_estimate method, this will its results. This can be used to convert the parameter vector to a more useful form, like a fatiando.mesher object.

fit()¶

Solve for the parameter vector that minimizes this objective function.

Uses the optimization method and parameters defined using the config method.

The estimated parameter vector can be accessed through the p_ attribute. A (possibly) formatted version (converted to a more manageable type) of the estimate can be accessed through the property estimate_.

fmt_estimate(p)¶

Called when accessing the property estimate_.

Use this to convert the parameter vector (p) to a more useful form, like a geometric object, etc.

Parameters:

• p

  : 1d-array

  The parameter vector.

Returns:

• formatted

  Pretty much anything you want.

gradient(p)¶

The gradient vector of the misfit function.

\[\bar{g} = -2\bar{\bar{J}}^T\bar{r}\]

where \(\bar{\bar{J}}\) is the Jacobian matrix and \(\bar{r}\) is the residual vector.

Parameters:

• p

  : 1d-array

  The parameter vector where the gradient is evaluated

Returns:

• gradient

  : 1d-array

  The gradient vector.

hessian(p)¶

The Hessian of the misfit function with respect to the parameters.

Calculated using the Gauss approximation:

\[\bar{\bar{H}} \approx 2\bar{\bar{J}}^T\bar{\bar{J}}\]

where \(\bar{\bar{J}}\) is the Jacobian matrix.

For linear problems, the Hessian matrix is cached in memory, so calling this method again will not trigger a re-calculation.

Parameters:

• p

  : 1d-array

  The parameter vector where the Hessian is evaluated

Returns:

• hessian

  : 2d-array

  The Hessian matrix

regul_param¶

The regularization parameter (scale factor) for the objetive function.

Defaults to 1.

residuals(p=None)¶

Calculate the residuals vector (observed - predicted data).

Parameters:

• p

  : 1d-array or None

  The parameter vector used to calculate the residuals. If None, will use the current estimate stored in estimate_.

Returns:

• residuals

  : 1d-array or list of 1d-arrays

  The residual vector. If this is the sum of 1 or more Misfit instances, will return the residual vector from each of the summed misfits in the order of the sum.

set_weights(weights)¶

Set the data weights.

Using weights for the data, the least-squares data-misfit function becomes:

\[\phi = \bar{r}^T \bar{\bar{W}}\bar{r}\]

Parameters:

• weights

  : 1d-array or 2d-array or None

  Weights for the data vector. If None, will remove any weights that have been set before. If it is a 2d-array, it will be interpreted as the weight matrix \(\bar{\bar{W}}\). If it is a 1d-array, it will be interpreted as the diagonal of the weight matrix (all off-diagonal elements will default to zero). The weight matrix can be a sparse array from scipy.sparse.

value(p)¶

Calculate the value of the misfit for a given parameter vector.

The value is given by:

\[\phi(\bar{p}) = \bar{r}^T\bar{\bar{W}}\bar{r}\]

where \(\bar{r}\) is the residual vector and \(bar{\bar{W}}\) are optional data weights.

Parameters:

• p

  : 1d-array or None

  The parameter vector.

Returns:

• value

  : float

  The value of the misfit function.

fatiando.geothermal.climsig.abrupt(amp, age, zp, diffus=31.5576)[source]¶

Calculate the residual temperature profile in depth due to an abrupt temperature perturbation.

Parameters:

• amp

  : float

  Amplitude of the perturbation (in C)

• age

  : float

  Time since the perturbation occured (in years)

• zp

  : array

  Arry with the depths of computation points along the well (in meters)

• diffus

  : float

  Thermal diffusivity of the medium (in m^2/year)

See the default values for the thermal diffusivity in fatiando.constants.

Returns

• temp

  : array

  The residual temperatures measured along the well

fatiando.geothermal.climsig.linear(amp, age, zp, diffus=31.5576)[source]¶

Calculate the residual temperature profile in depth due to a linear temperature perturbation.

Parameters:

• amp

  : float

  Amplitude of the perturbation (in C)

• age

  : float

  Time since the perturbation occured (in years)

• zp

  : array

  The depths of computation points along the well (in meters)

• diffus

  : float

  Thermal diffusivity of the medium (in m^2/year)

See the default values for the thermal diffusivity in fatiando.constants.

Returns

• temp

  : array

  The residual temperatures measured along the well