Method of calculating temperature of a geological structure

Abstract

A method of calculating the temperature of a geological structure is disclosed, wherein there is provided a magnetic parameter of the geological structure. The method includes inverting the magnetic parameter to estimate the temperature of the geological structure.

Claims

1. A method of prospecting for geothermal energy, petroleum, and/or minerals comprising calculating a temperature of a geological structure and using the temperature in a decision-making process for drilling a well or for mining a mine; wherein there is provided at least one magnetic parameter of the geological structure, the method comprising: inverting only the at least one magnetic parameter to estimate the temperature of the geological structure; wherein the at least one magnetic parameter comprises a total magnetization of the geological structure.

2. The method as claimed in claim 1, herein the total magnetization of the geological structure comprises both an induced magnetization and a remnant magnetization of the geological structure.

3. The method as claimed in claim 1, further comprising providing magnetic data and inverting the magnetic data to provide the at least one magnetic parameter.

4. The method as claimed in claim 3, wherein the magnetic data comprises magnetic potential field data.

5. The method as claimed in claim 3, wherein inverting the magnetic data to provide the at least one magnetic parameter is performed using a map inversion method.

6. The method as claimed in claim 3, wherein inverting the magnetic data to provide the at least one magnetic parameter provides laterally varying magnetization averaged over a depth interval.

7. The method as claimed in claim 3, wherein the magnetic data is acquired with a magnetometer.

8. The method as claimed in claim 1, wherein inverting the at least one magnetic parameter to estimate the temperature of the geological structure is performed using a Bayesian inversion method and a phenomenological model.

9. The method as claimed in claim 8, wherein the phenomenological model provides a relationship between the at least one magnetic parameter and the temperature of the geological structure.

10. The method as claimed in claim 8, wherein the phenomenological model uses or is based on: (i) a spinodal exsolution model and/or an Ising model from quantum statistics; or (ii) a bi-nodal exsolution model and/or JMAK kinetics.

11. The method as claimed in claim 8, wherein the phenomenological model comprises one or more parameters and the method comprises calibrating the one or more parameters.

12. The method as claimed in claim 11, wherein the one or more parameters comprise an initial titanium fraction of lavas at a time of deposition and/or a total percentage of magnetic material in a subsurface.

13. The method as claimed in claim 1, further comprising providing estimates of Curie temperatures, mineral composition, and/or fractions of different titanomagnetite phases.

14. The method as claimed in claim 1, further comprising selecting a phenomenological model that defines a relationship between the at least one magnetic parameter and the temperature of the geological structure.

15. The method as claimed in claim 1, further comprising finding an uncertainty in the calculated temperature.

16. The method as claimed in claim 1, wherein the temperature is calculated for a specific point/location/volume/space of the geological structure, the specific point/location/volume/space corresponding to a point/location/volume/space of the at least one magnetic parameter.

17. A method of producing a temperature model of the geological structure comprising performing the method of claim 1.

18. A method of harnessing the geothermal energy comprising performing the method of claim 1.

19. A method of obtaining the petroleum from the geological structure comprising performing the method of claim 1.

20. A method of mining for the minerals from the geological structure comprising performing the method of claim 1.

21. A computer program product comprising computer readable instructions that, when run on a computer, is configured to cause a processer to perform the method of claim 1.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) Preferred embodiments of the invention will now be discussed, by way of example only, with reference to the accompanying drawings, in which:

(2) FIG. 1 is the general multi-geophysical Bayesian network for estimation of the subsurface temperature T; and

(3) FIG. 2 is the Bayesian network for estimation of subsurface temperature T from magnetic anomaly data BA.

DETAILED DESCRIPTION

(4) An embodiment of the present invention provides a method of calculating or estimating the temperature of a geological structure. The method comprises two main steps: (i) inverting magnetic data to provide the magnetization of the geological structure; and (ii) inverting the magnetization to estimate the temperature of the geological structure.

(5) As discussed above, magnetic data has not previously been used to compute subsurface temperature.

(6) Embodiments of the present invention make use of the magnetic branch of the Bayesian network shown in FIG. 1 (and discussed above), including both induced and remnant magnetization.

(7) The magnetic branch of the Bayesian network is shown in more detail in FIG. 2. Here, T.sub.S is surface temperature, A.sub.TIR is the thermal infrared band from satellite data (LANDSAT-8), u.sub.0 and u are initial and present-day fractions of ulvöspinel, δ is the fraction of titanomagnetite oxidized to titanomaghemite, T.sub.C is Curie temperature, B.sub.0 is the earth magnetic background field (magnitude, inclination and declination), M is total magnetization (induced and remnant), B.sub.A is magnetic potential field data and c.sub.0 is the overall fraction of titanomagnetite in the basaltic rock.

(8) In principle, all of the variables in the Bayesian network of FIG. 2 can be regarded as stochastic. However, here, a simplified approach is used, taking only B.sub.A, M and T as stochastic variables. The other variables are treated as deterministic hyperparameters, or as having a delta-function distributions. Thus, the main parameters of interest here are magnetic data B.sub.A, rock magnetization M, and subsurface temperature T.

(9) The Bayesian network can be applied to obtain the joint distribution for a set of parameters, incorporating the principle of conditional independence. The joint probability of a set of stochastic nodes {x.sub.1, . . . , x.sub.n} can be written as

(10) $\begin{array}{cc}p\left(x_1,.Math.,x_n\right)=\underset{i}{.Math.}p\left(x_i.Math.x_i^{\mathrm{pa}}\right),& \left(1\right)\end{array}$
where x.sup.pa.sub.i denotes the parents of x.sub.i, i.e. nodes on the level above in the network.

(11) Using equation (1), and marginalizing hidden variables, the posterior distribution for subsurface temperature T given magnetic anomaly data BA can be written as
p(T|B.sub.A)=C∫dMp(T|M)p(M|B.sub.A)p(T)  (2)
where C is the normalization factor, and M is the magnitude of the magnetization. The integral marginalizes the magnetization M.

(12) As explained above, in practice, the inversion is performed in two separate steps: (i) the magnetic data BA are inverted to calculate the total magnetization M; then (ii) the total magnetization M is inverted to determine the subsurface temperature T.

(13) At step (i), magnetization M is computed by inversion of magnetic anomaly data B.sub.A. Using Bayes rule, the following is obtained
p(M|B.sub.A)∝p∫(B.sub.A|M)p(M),  (3)

(14) This step is linear if magnetic properties constrained on geometry are inverted for. This is discussed below.

(15) At step (ii), subsurface temperature T is computed by inversion of magnetization M. Again using Bayes rule, the following is obtained
p(T|M)∝p(M|T)p(T).  (4)

(16) This involves a non-linear phenomenological relationship between magnetization and temperature, which is discussed below.

(17) Finally, the posterior distribution p(T|B.sub.A) is obtained by means of equation (2). The marginalization of M can be written on a convolution form, which allows for fast and efficient computation using the fast Fourier transform (FFT).

(18) The magnetic data B.sub.A are magnetic potential field data measured with a magnetometer.

(19) Step (i) is performed using a map inversion method, e.g. a Marquardt-Levenberg type map inversion method, to determine laterally varying total magnetization averaged over a relevant depth interval. This is described in more detail below. Alternatively, step (i) can be performed using any other (e.g. 3D) inversion method serving the same purpose, such as the methods presented by Li and Oldenburg (1996) or Hokstad et al. (2017).

(20) The magnetic data B.sub.A can be acquired prior to any drilling and, in that case, can be used to provide a pre-drilling estimate of the temperature by performing the method according to the invention. Additionally or alternatively, magnetic data B.sub.A can be acquired during (or after) the drilling of a well. These data can be used to provide a during-drilling or post-drilling estimate of the temperature by performing the method according to the invention, which may be considered to be an update to the existing estimates.

(21) The map inversion method used in step (i) is now described in more detail.

(22) Magnetic anomaly data are due to dipole sources in the subsurface. Therefore, the anomalous magnetic field B can be obtained by double differentiation of an inverse distance potential

(23) $\begin{array}{cc}B\left(x\right)=-\frac{{\mu }_0}{4\pi }\nabla \nabla .Math.\int \frac{M\left(x^{\prime }\right)}{.Math.x^{\prime }-x.Math.}{\mathrm{dV}}^{\prime },& \left(5\right)\end{array}$
where M is the magnetization vector (including induced and remnant magnetization), x′ denotes the positions of magnetic sources, x is the position where the magnetic field is measured, and μ.sub.0 is the vacuum permeability.

(24) In standard acquisition of magnetic potential field data, only the magnitude of the total magnetic field, the total magnetic intensity is obtained
B.sub.TMI=|B.sub.0+B|=√{square root over (|B.sub.0|.sup.2+|B|.sup.2+2B.sub.0.Math.B)},  (6)
where B.sub.0 is the earth magnetic background field.

(25) The total magnetic anomaly B.sub.A is obtained by subtracting the local magnetic field,
B.sub.A=|B.sub.0+B|−|B.sub.0|,  (7)

(26) Expanding the square-root in equation (7) to second order in B, the total magnetic anomaly can be written as

(27) $\begin{array}{cc}{B}_A\simeq .Math.B.Math.\cos \theta +\frac{{.Math.B.Math.}^2{\sin }^2\theta }{.Math.B_0.Math.},& \left(8\right)\end{array}$
where θ is the angle between B and B.sub.0.

(28) Introducing the unit vector t=B.sub.0/|B.sub.0|, in the direction of B.sub.0, the TMA can be written as

(29) $\begin{array}{cc}{B}_A\simeq t.Math.B+\frac{{.Math.t\times B.Math.}^2}{.Math.B_0.Math.}.& \left(9\right)\end{array}$
When |B|<<|B.sub.0|, the second term is insignificant, and the total magnetic anomaly is essentially the projection of the anomalous field onto the direction of the earth magnetic background field.

(30) This assumption is used in the following.

(31) Magnetic data capture lateral changes in magnetization very well, but have very little depth resolution. Therefore, 3D magnetic inversion is severely ill-posed, and must be heavily constrained to give meaningful results. The depth distribution of magnetic anomalies from 3D inversion is to a large degree controlled by the regularization used, and should not be over-interpreted.

(32) In the present invention, this problem is addressed using a map-inversion method, as described in the following.

(33) Magnetization is assumed to be independent of depth, or it represents an average over the depth range z.sub.1 to z.sub.2, such that

(34) $\begin{array}{cc}M\left(x^{\prime },y^{\prime }\right)=\frac{1}{z_2-z_1}{\int }_{z_1}^{z_2}M\left(x^{\prime },y^{\prime },z^{\prime }\right)& \left(10\right)\end{array}$
Then the integral over z′ in equation 5 can be carried out, giving

(35) $\begin{array}{cc}{{\int }_{z_1}^{z_2}\frac{1}{r}{\mathrm{dz}}^{\prime }=\ln \left[\left(z^{\prime }-z\right)+r\right].Math.}_{z_1}^{z_2}=F\left(z_2\right)-F\left(z_1\right).& \left(11\right)\end{array}$
Where, for later convenience, we have introduced the distance
r=√{square root over ((α′−x).sup.2+(y′−y).sup.2+(z′−z).sup.2)},  (12)
and the function
F(z′)=ln[(z′−z)+r].  (13)
Substituting in equation 5, and using equation 9, the total magnetic anomaly can be written as

(36) $\begin{array}{cc}{B}_A\left(x\right)\simeq t.Math.B\left(x\right)=\frac{{\mu }_0}{4\pi }\underset{i=1}{\overset{3}{.Math.}}\underset{j=1}{\overset{3}{.Math.}}\int t_i{A}_{\mathrm{ij}}\left(x,x^{\prime }\right)M_j\left(x^{\prime },y^{\prime }\right){\mathrm{dx}}^{\prime }{\mathrm{dy}}^{\prime },& \left(14\right)\end{array}$
where t.sub.i and M.sub.j are elements of the vectors t and M, respectively, and A.sub.ij are the elements of a second-order tensor given by

(37) $\begin{array}{cc}{A}_{\mathrm{ij}}=\frac{{\partial }^{2}F\left(z_2\right)}{\partial x_i\partial x_j}-\frac{{\partial }^{2}F\left(z_1\right)}{\partial x_i\partial x_j},& \left(15\right)\end{array}$
where the derivatives are with respect to horizontal coordinates (x.sub.1, x.sub.2)=(x, y), and

(38) $\begin{array}{cc}\begin{array}{c}\frac{{\partial }^{2}F\left(z^{\prime }\right)}{\partial x_i\partial x_j}=-\frac{1}{{.Math.\left(z^{\prime }-z\right)+r.Math.}^2}\left[{\delta }_{i3}+\frac{\left(x_i^{\prime }-x_i\right)}{r}\right]\left[{\delta }_{j3}+\frac{\left(x_j^{\prime }-x_j\right)}{r}\right],\\ +\frac{1}{.Math.\left(z^{\prime }-z\right)+r.Math.}\left[\frac{{\delta }_{\mathrm{ij}}}{r}-\frac{\left(x_i^{\prime }-x_i\right)\left(x_j^{\prime }-x_j\right)}{r^3}\right],\end{array}& \left(16\right)\end{array}$
and δ.sub.ij is the Kronecker delta function.

(39) From the forward model equation 14, it is clear that three-component magnetization M cannot be readily obtained by inversion of scalar data B.sub.A. However, if an estimate of the direction of magnetization is known, it may be attempted to invert for the magnitude, which is a feasible problem.

(40) Hence, it can be written that
M=Me,  (17)
where M(x′, y′) and e are the magnitude and a unit vector in the direction of the magnetization, respectively. In many cases the magnetization is approximately parallel to the paleomagnetic field at the time when the rock solidified. For very young igneous rocks it is reasonable to assume e=t. When the magnetization M is parallel to the inducing field B.sub.0, it is difficult to distinguish remnant magnetization M.sub.R from induced magnetization M.sub.I=χH.sub.0, where H.sub.0=B.sub.0/μ.sub.0. Without additional information, the two parts are often mixed together in an apparent susceptibility
χ.sub.eff=M/H.sub.0=χ+M.sub.R/H.sub.0.  (18)

(41) In general, a tectonic reconstruction and kinematic restoration may be needed to estimate e.

(42) Assuming a geomagnetic problem with n.sub.d measurements and n.sub.m model parameters, the forward problem can be written on matrix-vector form as
d=Lm+n,  (19)
where d is a n.sub.d×1 data vector with component d.sub.i=B.sub.A(x.sub.i,y.sub.i) is a n.sub.m×1 model vector with components m.sub.j=M(x.sub.j,y.sub.j), n is noise, and L is a n.sub.d×n.sub.m matrix. Assuming that e=t, the elements of the matrix L become

(43) 0$\begin{array}{cc}{L}_{\mathrm{ij}}=\frac{{\mu }_0}{4\pi }\underset{k=1}{\overset{3}{.Math.}}\underset{l=1}{\overset{3}{.Math.}}{t}_k{A}_{\mathrm{kl}}\left(x_i,x_j\right)t_l.& \left(20\right)\end{array}$

(44) In a statistical formulation of the inverse problem, using Bayes rule, the posterior distribution p(m|d)| can be written as
p(m|d)=Cp(d|m)p(m)=Ce.sup.−Φ,  (21)
where C is a normalization factor, and Φ is the objective function for an optimization-type inversion,
Φ=½|d−Lm|.sup.TΣ.sub.d.sup.−1|d−Lm|+½[m].sup.TΣ.sub.m.sup.−1m.sup.Tm.  (22)

(45) It is assumed that the prior distribution for m is Gaussian with zero mean, which is reasonable for inversion of anomaly data. Since the forward model equation 19 is linear in m, the likelihood and posterior distributions are Gaussian. Then maximizing the posterior probability is equivalent to minimizing the objective function Φ. This can be achieved with the least squares-pseudo inverse,
μ.sub.m|d≅m=|L.sup.TΣ.sub.d.sup.−1L+Σ.sub.m.sup.−1|.sup.−1L.sup.TΣ.sub.d.sup.−1d,  (23)
where μ.sub.m|d denotes the posterior mean of the magnetization M. For simplicity, Σ.sub.d is chosen to be proportional to the identity matrix. If Σ.sub.m.sup.−1 is chosen as
Σ.sub.m.sup.−1=λdiag(L.sup.TL),  (24)
the classical Marquardt-Levenberg solution is obtained. If Σ.sub.m.sup.−1 is chosen as
Σ.sub.m.sup.−1=λD.sup.TD,  (25)
where D is the first-derivative matrix, a smoother solution is obtained. Here, λ is a parameter controlling the relative weights put on the data term and regularization term of the objective function.

(46) The posterior covariance matrix is independent of the data, and can be written explicitly as
Σ.sub.m|d=Σ.sub.m−Σ.sub.mL.sup.T[LΣ.sub.mL.sup.T+Σ.sub.d].sup.−1LΣ.sub.m.  (26)
The posterior variances are the diagonal elements of Σ.sub.m|d. The second term on the right hand side is positive definite. Hence, the effect of adding information from data is to reduce the uncertainty.

(47) Step (ii) uses a Bayesian inversion method involving a phenomenological model relating the magnetization to the temperature, to calculate the temperature value from the magnetization. The phenomenological model is based on spinodal exsolution of titanomagnetite and the 1D Ising model from quantum statistics. This is described in more detail below. As a by-product of this step, estimates of the Curie temperatures and fractions of different titanomagnetite phases are also obtained.

(48) In order to obtain a spatially dependent 3D temperature function, the temperature of the geological structure is calculated point-wise for multiple different points/locations/volumes/spaces in the geological structure. As can be appreciated, the magnetic parameter(s) may vary over the space of the geological structure, and this may correspond to a spatially varying temperature.

(49) The phenomenological model relating magnetization to underlying rock properties is now described.

(50) The magnetic properties of rocks are mainly determined by a small number of strongly magnetic minerals. Most important are titanomagnetite, which may subsequently oxidize to titanomaghemite. Still, titanomagnetite only makes up a relatively small fraction of the total rock volume, typically 1-10%. The remnant and induced magnetization depend strongly on temperature, in addition to chemical processes and hydrothermal alteration. Magnetization is strongly temperature dependent, and can be pictured as a battle between quantum-mechanical order and thermal disorder. When magma cools and solidifies, single-domain titanomagnetite forms as a solid solution of magnetite and ulvöspinel. Upon further cooling, this solid solution may enter the miscibility gap, where it becomes unstable, and decomposes into new phases with different composition.

(51) To describe the thermal effects on magnetization, a simple model is needed that statistically links quantum effects to the macroscopic magnetization. Here, a simplified version of the Ising (1925) model from quantum statistics is used. In the Ising model, the Hamiltonian (energy) for a system of particles with interacting spins σ.sub.j, and under the influence of an external magnetic field H is given by

(52) $\begin{array}{cc}\mathscr{H}=-J\underset{.Math.\mathrm{ij}.Math.}{.Math.}{\sigma }_i{\sigma }_j-{\mu }_{H}H\underset{j}{.Math.}{\sigma }_j,& \left(27\right)\end{array}$
where μ.sub.B is the magnetic moment of each particle, J is the inter-particle interaction strength, and the notation <ij> means that the sum runs over all nearest-neighbors combinations i; j. In equation 27, the first term represents remnant magnetization, and the second term induced magnetization.

(53) After some algebra, using methods of statistical mechanics and the Boltzmann distribution, the macroscopic magnetization of an ensemble of particles can be written as

(54) $\begin{array}{cc}\frac{M}{M_0}=\tanh \left[\frac{T_{C}M}{{\mathrm{TM}}_0}+\frac{\mathrm{CH}}{{\mathrm{TM}}_0}\right].& \left(28\right)\end{array}$
Here M is the magnitude of the magnetization, M.sub.0 is the magnetization at absolute zero, T is the temperature (in Kelvin) T.sub.C is the Curie temperature, and C is the Curie constant. The derivation of equation 28 is outlined below. Equation 28 is a transcendental equation which can be solved graphically or numerically. The susceptibility χ is given by

(55) $\begin{array}{cc}\chi =\frac{M\left(T,H\right)-M\left(T,0\right)}{H},& \left(29\right)\end{array}$
where χ is the magnetic susceptibility, and M(T, 0)=M.sub.R is the remnant magnetization. For large temperatures, equation 28 can be linearized, which leads to the Curie-Weiss law,

(56) $\begin{array}{cc}\chi =\frac{M}{H}=\frac{C}{T-T_C}.& \left(30\right)\end{array}$

(57) The Ising model can be used to represent the magnetic properties of a homogeneous single-domain magnetic phase. For a heterogeneous mix of magnetic phases, a linear combination of single-phase models is used,

(58) $\begin{array}{cc}M\left(T,H\right)=\underset{j}{.Math.}{w}_j{M}_j\left(T,H\right),& \left(31\right)\end{array}$
where w.sub.j is the weight of magnetic phase j. The linear superposition is supported by the linearity in M of the equations 5 and 15.

(59) With equation 28, only positive magnetization can be modelled. Negative remnant magnetization is, however, frequently observed in magnetic data, also for young MORB. Negative magnetization can be modelled by means of self-reversal for one or more of the magnetic phases. In practice this is achieved by means of negative weight w.sub.i, and reversing the sign of H in equation 28.

(60) The model given in equation 31 can represent a large variety of temperature dependence of magnetization, by superposition of two or three mechanisms with different Curie temperatures, similar to those shown in experimental studies.

(61) When an igneous rock forms by cooling of magma, titanomagnetite forms as a solid solution of ulvöspinel and magnetite. Also, a part of the titanomagnetite may be transformed by high temperature oxidation to titanomaghemite.

(62) The spinodal decomposition takes place for a simple two-component solution. A homogeneous solid solution is formed at temperature T.sub.0 with initial fractions u.sub.0 and v.sub.0 of the two components. Upon subsequent cooling, to temperature T the solid solution enters the miscibility gap, where it becomes unstable, and decompose into new phases with compositions u.sub.1, v.sub.1 and u.sub.2, v.sub.2. The phases may continue to decompose by bi-nodal decomposition to phases with compositions u.sub.3, v.sub.3 and u.sub.4, v.sub.4. However, this relies on the process of nucleation, and is a much slower process. From mass conservation, u and v are not independent, and the system can be described by one of them. In the following, u is used.

(63) The spinodal decomposition is controlled by the Gibbs free energy. For the simplest model, the Gibbs free energy and spinodal are symmetric around u=½. Laboratory studies on spinodal decomposition of titanomagnetites show that the spinodal does not have this symmetry. Gibbs free energy for a two-component solution consists of an enthalpy term and an entropy term. In an external magnetic field, the entropy of the magnetite will be lower than the entropy of ulvöspinel, because magnetite is magnetic, whereas ulvöspinel is not. This asymmetry can be captured if the Gibbs free energy G of the system is expressed as
G(u,v,T)=W(uv+uδ+vδ)+k.sub.BT[(1+h)uln u+(1−h)vln v+δ ln δ]|,  (32)
where u is the fraction of magnetite, v is the fraction of ulvöspinel, and δ is the fraction of titanomaghemite. δ is assumed to be small, and constant during the spinodal decomposition. By mass conservation the conditions u+v+δ=1 and θu/θG=−1 are obtained. Moreover, W is the mixing enthalpy parameter, k.sub.B is the Boltzmann constant, and the parameter h accounts for the different magnetic entropies of magnetite and ulvöspinel in the presence of an external magnetic field.

(64) The first and second derivatives of the Gibbs free energy are obtained as

(65) $\begin{array}{cc}\frac{\partial G}{\partial u}=-W\left(u-v\right)+k_{B}T\left[\left(1+h\right)\ln u-\left(1-h\right)\ln v+2h\right],& \left(33\right)\end{array}$$\mathrm{and}$$\begin{array}{cc}\frac{{\partial }^{2}G}{\partial u^2}=-2W+k_{B}T\left[\frac{\left(1+h\right)}{u}+\frac{\left(1-h\right)}{v}\right]& \left(34\right)\end{array}$

(66) The spinodal decomposition takes place when

(67) $\frac{{\partial }^{2}G}{\partial u^2}<0,$
and terminates at the spinode defined by

(68) $\frac{{\partial }^{2}G}{\partial u^2}=0.$
Equation 34 is a square equation for u, which gives two solutions u.sub.1 and u.sub.2, for two different solid solutions with different magnetite and ulvöspinel fractions. Setting h=0 and δ=0, we obtain spinodal symmetric around u=½. The parameters W and h were adjusted to fit the spinodal model to experimental measurements. Given temperature, the model then predicts the equilibrium fractions of ulvöspinel and magnetite for the phases with composition u.sub.1 and u.sub.2, respectively. Moreover, from the empirical model presented by Lattard et al. (2006),
T.sub.C(u)=−150u.sup.2−580u+851,  (35)
the corresponding Curie temperatures of the two phases are also obtained.

(69) The main constituents of MORB are olivine, pyroxene and plagiochlase. Though titanomagnetite and sometimes titanomaghemite are the main minerals responsible for the magnetization of the basalt, these minerals will normally make up only a small percentage of the total rock volume, typically 2-10%.

(70) By stoichiometric balancing, the weight functions for the Ising model equation 31 are obtained as

(71) $\begin{array}{cc}{w}_1=c_0\left(1-\delta \right)\frac{u_0-u_2}{u_1-u_2}& \left(36\right)\end{array}$$\begin{array}{cc}{w}_2=c_0\left(1-\delta \right)\frac{u_1-u_0}{u_1-u_2}& \left(37\right)\end{array}$$\begin{array}{cc}{w}_3=c_0\delta & \left(38\right)\end{array}$
where c.sub.0 is the overall volume fraction of magnetic minerals (titanomagnetite and titanomaghemite) in the rock.

(72) Due to the low depth resolution of magnetic potential field data, a simplified representation in terms of depth-average magnetization was chosen, as given in equation 10. Consequently, a simplified parameterization for subsurface temperature vs depth was also sought. The subsurface temperature T is represented as
T(z,β)=T.sub.∞+(T.sub.S−T.sub.∞)e.sup.−βz  (39)
where T.sub.S is the land surface temperature, and T.sub.∞ is the asymptotic value for z.fwdarw.∞. With the decay parameter β varying from 0.1 to 0.5 km.sup.−1, anything from conductive thermal gradients of order 40° C./km to geothermal anomalies close to the boiling curve can be represented.

(73) The land surface temperature T.sub.S may vary significantly in active geothermal areas, e.g. near hot springs and fumaroles. A simple option is to express T.sub.S as a function of β
T.sub.S(β)=αβ.sup.n+b  (40)
where n=4 or n=6 give reasonable results.

(74) The land surface temperature Ts can also be constrained on observed data. The land surface temperature can be estimated using remote-sensing satellite data. Here, a method utilizing LANDSAT-8 bands 4,5 and 8, i.e. the red, near infrared (NIR) and thermal infrared (TIRS) bands is applied. The method of Avdan and Jovanovska (2016) has been adapted to a statistical inversion framework, similar to that presented for magnetization above. However, the resolution of the TIRS bands (90 m) is not able to capture the high temperatures near hot springs.

(75) The decay parameter β can be associated with an average temperature T(β) defined by

(76) 0$\begin{array}{cc}\hat{T}\left(\beta \right)=\frac{1}{z_2-z_1}{\int }_{z_1}^{z_2}T\left(z,\beta \right)\mathrm{dz}& \left(41\right)\end{array}$

(77) Given average magnetization from inversion of magnetic anomaly data, the average magnetization can be inverted for average subsurface temperature of the magnetic source layer. Using equations 39 and 41, the link between average temperature and temperature vs depth is established at the forward modelling stage of the Bayesian inversion, i.e. as part of computing the likelihood part of the posterior probability distribution.

(78) The above method can be used when prospecting for geothermal energy, petroleum, and/or minerals (e.g. metal sufides), e.g. when planning and performing geothermal/petroleum well drilling or mineral mining operations.

(79) In one embodiment, the calculated temperature is used prior to drilling or mining, when deciding where to drill the well and/or how deep to drill the well, or where/how deep to mine the mine.

(80) In the same or other embodiments, the calculated temperature is used during or after the drilling of the well or mining of the mine, when deciding in which direction or to which depth to continue drilling or mining.

(81) The temperature estimate can be updated during drilling/mining based on new measure data.

(82) It should be apparent that the foregoing relates only to the preferred embodiments of the present application and the resultant patent. Numerous changes and modification may be made herein by one of ordinary skill in the art without departing from the general spirit and scope of the invention as defined by the following claims and the equivalents thereof.

REFERENCES

(83) Advan, U., and G. Jovanovska, 2016, Algorithm for automated mapping pf land surface temperature using LANDSAT 8 satellite data; Journal of Sensors, 2016, 1-8. Hokstad, K., and K. Tänaysuu-Mileviciene, 2017, Temperature prediction by multigeophysical inversion: Application to the IDDP-2 well at Reykjanes, Iceland: GRC Transactions, 42, 1141-1152. Lattard, D., R. Engelmann, A. Kontny, and U. Sauerzapf, 2006, Curie temperatures of synthetic titanomagnetites in the Fe—Ti—O system: Effects of composition, crystal chemistry, and thermomagnetic methods: J. Geophys. Res., 111. Li, Y., and D. Oldenburg, 1996, 3-D inversion of magnetic data: Geophysics, 61, 394-408.