Radiogenic heat production in the crust from inversion of gravity and magnetic data

Understanding heat flow and its lateral variability is important for petroleum prospecting. We present a methodology and workflow for estimation of radiogenic heat production in the crust and heat flow from geophysical data. The surface heat flow is determined by three main components: heat from Earth’s interior, heat produced in the crust, and heat produced in the sediments. To estimate the crustal heat contribution, we have developed a statistical inversion approach. We first estimate density and susceptibility by inversion of crustal gravity and magnetic anomalies. Then we compute the radiogenic heat production in the crystalline crust by Bayesian rock-physics inversion of density and susceptibility. We demonstrate the proposed methodology with two case examples. In the first example, we compute radiogenic heat production from density and susceptibility measured on rock samples from onshore Norway. The estimated average value of the radiogenic heat production in the rock samples is about 2 μW/m3, with uncertainty of ± 0.5 μW/m3. The second example is from the greater Barents Sea, where we applied a full workflow including a large-scale lithospheric modelling to estimate the mantle heat flow, 3D gravity and magnetic inversion, and rock-physics inversion. The average heat production in the Barents Sea crust is about 1.5 μW/m3. The predicted heat production and heat flow are in good agreement with results presented by other authors.


Introduction
The thermal state and thermal history of the subsurface are of major importance for the petroleum prospectivity of sedimentary basins.The thermal history controls both the maturation of source rocks and the quality of reservoir rocks.In basin modelling, temperature and heat-flow data are needed for calibrating and constraining models.Standard calibration data are bottom-hole temperatures and vitrinite reflectance, which are only available in boreholes.In frontier areas, boreholes are often scarce.In some cases, diagenetic surfaces like opal A-CT, can be used as thermal indicators.Heat-flow measurements from piston coring may be available, but are often as uncertain as bottomhole temperatures (Khuturskoi et al., 2008).
In addition to thermal conductivity of the lithospheric layers, the temperature distribution in a sedimentary basin is determined by three major components: Radiogenic heat production is caused by decay of the long-lived radiogenic isotopes of thorium, uranium and potassium (e.g., Turcotte & Schubert, 2002), mainly in the continental crust.Oceanic crust and lithospheric mantle do not contain significant amounts of radioactive isotopes (Hasterok, 2010;Allen & Allen, 2013).Therefore, it is essential to find the RHP sources within the continental crust and sediments.The RHP in the sediments can be constrained by well data (preferably spectral gamma-ray logs), while to estimate the RHP in the crust on regional scales, various rocks-physics parameters and data need to be utilised as a proxy for direct measurements.Hence, we have developed a method using geophysical data as an information carrier between wells to constrain present-day heat flow.Rybach (1978) presented an exponential relationship between RHP and seismic velocity.Rybach & Buntebarth (1984) demonstrated a correlation for both seismic velocity and density with radiogenic heat production.Adding to this, we have developed a rock-physics model correlating RHP to magnetic susceptibility.The correlation between heat production and geophysical parameters has been disputed.Kukkonen & Peltoniemi (1998) investigated 2700 core samples from Finland, and concluded that the correlation between heat production and other geophysical and petrophysical parameters were generally low.They argued that heat production is due to the presence of thorium, uranium and potassium in the rocks, whereas geophysical parameters such as density, velocity, susceptibility and resistivity are controlled by other factors.
On the other hand, a correlation between RHP and geophysical parameters was provided by Mernagh & Miezitis (2008) and Hasterok & Webb (2017), based on geochemical analysis.Continental crust is formed through several stages of refinement, involving partial melting and fractionation of melt.In this process, thorium and uranium isotopes tend to follow the quartz, which is enriched in the upper part of the crust.Thorium and uranium are also sometimes referred to as 'incompatible elements' .They are passed on to the part of the melt with the lowest melting point (Mernagh & Miezitis, 2008).The melting temperature of the rock decreases with increasing quartz (SiO 2 ) content.
In previous work (Hokstad, 2014;Hokstad et al., 2014;Duffaut et al., 2017), we addressed the thermal conductivity of the sedimentary sequence, and the relationship between the thermal conductivity and seismic velocity.Here, we investigate the basal and surface heat flow.Our main objective is to use geophysical data and 3D geophysical models to constrain heat production and heat flow (Hokstad, et al., 2016).
In the simplest 1D steady-state approximation, the temperature distribution of the subsurface is given by Fourier's law, dT = q dz k ' (1) where dT/dz is the thermal gradient, q is the heat flow, and k is the thermal conductivity.The thermal gradient is invariant under the scaling q → cq and k → ck, where c is an arbitrary constant.Hence, when calibrating and comparing thermal modelling results, both heat flow and thermal conductivity should be considered.
We compute the mantle contribution to the heat flow using the LitMod3D modelling tool (Afonso et al., 2008, Fullea et al., 2009).To estimate radiogenic heat production in the crust, and its contribution to heat flow, we apply a novel multi-geophysical inversion method.First, we perform single-domain inversions of gravity and magnetic data.Second, we use a Bayesian inversion method to estimate RHP from 3D susceptibility and density models.Third, we integrate the RHP over the thickness of the crust, and add it to the mantle contribution.Finally, an estimate of surface heat flow can be obtained by adding heat production in the sedimentary package, e.g., based on gamma-ray logs (Rybach, 1986;Bücker & Rybach, 1996).
We demonstrate the proposed method in two case studies, onshore and offshore, respectively.In the first step we work with samples and use them to adjust the rock-physics model, as well as to test the inversion on the onshore samples from Norway.Then we apply the full workflow to gravity and magnetic data from the greater Barents Sea.The estimated basal and surface heat flow can be applied directly to compute temperature within the steady-state approximation.
Our approach differs from other published forward modelling approaches, such as Guillou et al. (1994), Zeyen & Fernàndez (1994), Afonso et al. (2008) and Fullea et al. (2009), by explicitly inverting for the RHP in the crust based on a specific work flow, rather than by assigning a certain value of the RHP to the crustal layers to match the measured surface heat flow.

Methods and general workflow
To obtain realistic surface heat flow, with lateral variations, all three main components mentioned above must be separately assessed.The RHP in the sediments is mainly constrained by gamma-ray logs and our in-house data confirm an average value of ~1 μW/m 3 .The heat produced within the crystalline crust (i.e., the crustal contribution) and the mantle heat flow estimated at the Moho (i.e., base of crust) are discussed in detail below.
the mantle can be modelled, since LitMod3D considers all the aforementioned geophysical parameters and data consistently.The modelled densities are thoroughly determined even in tectonic settings with pronounced undulations of the LAB.If an optimal solution cannot be achieved due to limited time and/or lack of constraining data, it is at least possible to test various scenarios of published depth to the LAB and exclude models that are not compatible with measured data.For this purpose, not only the global satellite gravity data are suitable, but also digital elevation models (DEM) have been shown to be extremely useful.The calculated elevation is strongly dependent on the lithospheric mantle structure and tectonic evolution of the area.

Heat production in the crust
We developed a novel method to compute the presentday heat production in the crystalline crust from susceptibility and density models obtained by inversion of potential field data.We assume that the relationship between geophysical data and RHP can be modelled in terms of a Bayesian network, or directed acyclic graph (DAG), as shown in Fig. 1.In this network, different geophysical data and models are conditionally independent, but depend on a common parameter, A, which for our problem is the RHP.From the Bayesian network, a statistical inversion scheme can be derived.
The inversion is performed in two steps.First, we compute geophysical models by inversion of geophysical

Mantle heat flow
The mantle heat flow in stable continental regions (shields and cratons) varies globally from 11 to 25 mW/ m 2 (Mareschal & Jaupart, 2013).The value of 18 mW/ m 2 has been found to be the best-fitting average Moho heat flow for these regions (Jaupart & Mareschal, 2007).However, the study by Hamza & Vieira (2012) shows a much higher mantle heat flow, varying from 20-40 mW/ m 2 on continents to more than 40 mW/m 2 in offshore areas (fig.7 in Hamza & Vieira, 2012).According to Turcotte & Schubert (2002), the average continental mantle heat flow is 28 mW/m 2 , and according to Allen & Allen (2013) it is 30 mW/m 2 .The basins of our interest are commonly located in the transitional onshoreoffshore regions, which are characterised by pronounced differences in the structures and are the least 'stable' or predictable based on global averages only.
Given the wide range of possible mantle heat flow, it is clear that more accurate values must be determined for each specific region of interest.For this purpose, we apply the LitMod3D modelling tool (Afonso et al., 2008;Fullea et al., 2009), which makes simultaneous use of geophysical and petrophysical data.This is an academic software package, not previously applied to petroleum exploration.LitMod3D is a forward modelling approach enabling development of large-scale deep models, down to 400 km depth.In our workflow, the present-day lithospheric structure down to the lithosphere-asthenosphere boundary (LAB) must be defined.In LitMod3D, the LAB is a thermal boundary layer represented by a certain isotherm (usually between 1200°C and 1315°C).Other input parameters needed are the geometries of the crustal and mantle layers, their thermal properties (thermal conductivity and heat production), densities in the crust and composition of the underlying mantle part.
Based on the given structure and the input parameters, first the large-scale thermal distribution is calculated.
Using the given chemical composition of the mantle, densities and pressure are calculated iteratively.Densities in the mantle part are calculated in situ as a function of pressure and temperature by means of a coupled petrological modelling procedure.Once the density and temperature distributions are known, a number of geophysical quantities can be computed: surface heat flow, elevation, geoid, Bouguer and free-air gravity anomalies.These calculated responses are compared to the measured data and used to adjust the model.This procedure is repeated until a good fit to all the observed data is achieved.From a final model, only the contribution of the mantle part is extracted, since the estimates of the crustal contribution are made with greater detail in a separate approach, described below.
An important feature of the LitMod3D is the modelling of several geophysical datasets simultaneously.Therefore, both thermal and compositional density variations in data.Second, we estimate radiogenic heat production by inversion of geophysical model parameters.This step is a rock-physics inversion, implemented as a statistical scheme.
In practice, we take a pragmatic approach, performing a stepwise multi-geophysical inversion by combining separate single-domain inversions.In doing so, we can utilise geophysical models computed with different methods by various software providers, and produce results with fast turnaround.Density and susceptibility models from 3D gravity and magnetic inversion are the most important input models for regional studies.
Gravity and magnetic data are commonly available with 3D coverage over large areas, and play an important role in early-phase exploration.Seismic velocity cubes can also be utilised if available.
Gravity and magnetic data contain important information about crustal composition, but have hardly any depth resolution (Blakely, 1996).Consequently, the inversion of gravity and magnetic data is very ill-posed.
To overcome this problem, we have developed a purposemade 3D inversion scheme, where we invert for crustal properties given geometrical constraints, i.e., top of the crystalline crust (crystalline basement) and Moho.The details of the gravity and magnetic inversions are given in the Appendix.
Prior to the inversion, we need to isolate the crustal anomalies.For the gravity data, this means stripping off the effects of the water layer and sedimentary rocks, and correcting for lateral undulations in the Moho depth.
For the magnetic field, only a reduction to pole filter is applied.Ideally, we should also remove the remanent part of the magnetisation, but in practice this is often difficult.
The outputs from the geophysical inversions are input to a statistical rock physics inversion.From the Bayesian network in Fig. 1, the posterior distribution for radiogenic heat production A, given one or more geophysical model parameters m i , can be written as (Cowell et al., 2007;Lilleborge et al., 2015) (2 where C is the normalisation factor, and p(A) is the prior distribution for the radiogenic heat production.Different lithology-dependent priors can be used if we have relevant knowledge of the nature of the crust (oceanic or continental).For the problem considered here, m i = {r, χ}, where r is density and χ is magnetic susceptibility.Assuming Gaussian errors, each of the likelihood functions p(m i |A) in the product on the righthand side of Equation ( 2), can be written as where σ ei is the error variance for model parameter m i , and F i (A) are rock-physics models relating the heat productions to geophysical model parameters.Note that the rock-physics models are forward models.One example, the relation between heat production and seismic velocity presented by Rybach (1978), is shown in Fig. 2. Similar models can be obtained for heat production vs. density and susceptibility.We define the error variance as the RMS error between the model curve and calibration data (Fig. 2).One or more geophysical parameters can be used to compute the posterior distribution for radiogenic heat production.We can for instance use density alone.However, the product of two (or more) likelihood functions makes the posterior distribution narrower, which implies better posterior mean and smaller variance.The error variances can also be used to adjust the relative weighting of the density and susceptibility in the inversion, e.g., increasing the susceptibility variance σ eχ if significant remanence is suspected in the magnetic data.
Given the posterior distribution, the mean and variance of the radiogenic heat production can be computed from the definitions.The heat production from the rockphysics inversion is then integrated from the top of the crystalline crust (top basement) to the top of the lower crust to obtain the crustal heat-flow contribution.If we lack the top of the lower crust horizon (i.e., the upper boundary of the lower crust), we assume that the upper and middle crust constitutes 2/3 of the total crustal thickness (Beardsmore & Cull, 2001).Most of the rocks containing the radiogenic isotopes are in the upper to middle crust (e.g., Allen & Allen, 2013).Therefore, the integration is performed only over this interval.The line crust (crystalline basement) and Moho.The details of the gravity and magnetic e given in the Appendix.
nversion, we need to isolate the crustal anomalies.For the gravity data, this means the effects of the water layer and sedimentary rocks, and correcting for lateral n the Moho depth.For the magnetic field, only a reduction to pole filter is applied.
hould also remove the remanent part of the magnetisation, but in practice this is t.
from the geophysical inversions are input to a statistical rock physics inversion.
yesian network in Feil! Fant ikke referansekilden., the posterior distribution for The outputs from the geophysical inversions are input to a statistical rock physics inve From the Bayesian network in Feil! Fant ikke referansekilden., the posterior distributio radiogenic heat production , given one or more geophysical model parameters  / , c written as (Cowell et al., 2007;Lilleborge et al., 2015)  | 2 , … ,  4 =    / |   , were measured, from which the heat production can be calculated using the standard equation published by Rybach (1988) and Hasterok & Webb (2017).We will refer to this as observed RHP data (Fig. 3A).
Measured densities and susceptibilities were inverted, for comparison with observed RHP values (Fig. 3B).We used these data to calibrate the crustal rock-physics model of the inversion to match average observed heat-production values per lithology (Fig. 4), particularly for the common crustal rocks (i.e., granites, rhyolites, gabbros, granulites, eclogites).Some of the less common compositions, crustal contribution is then combined with the mantle part from LitMod3D, to obtain the basal heat flow.In order to estimate the surface heat flow, heat production in the sediments (~1 μW/m 3 ) must be added to the basal heat flow.

Heat production onshore Norway
The Geological Survey of Norway (NGU) has measured density and susceptibility on rock samples from many locations onshore Norway (Olesen et al., 2010).In addition, for approximately 450 samples provided by NGU, also uranium, thorium and potassium contents  such as mangerite, tonalite, anorthosite and migmatite, do not follow the rock-physics model of the inversion.Therefore, the inverted results have larger deviations for these samples.
Subsequently, the rock-physics inversion was applied to approximately 18,000 samples, classified as igneous and metamorphic rocks, ranging from felsic to ultramafic in composition (Fig. 5).Most of the samples have an estimated RHP of ~2 µW/m 3 or less.A few samples have values up to 6 µW/m 3 .The uncertainty (posterior variance) is in the order of ± 0.5 µW/m 3 .Some areas, for instance the Oslo Graben and Østfold (Iddefjorden granite), stand out with relatively high values.The inversion results obtained from the onshore samples compare well with the heat-production map published by Slagstad (2008).The heat production in the Lofoten Islands, however, is slightly overestimated.
Nevertheless, the main objective of the Bayesian multigeophysical inversion is its application to 3D density and susceptibility models, obtained from gravity and magnetic data.From measurements on rock samples to gravity and magnetic inversion grids is a big variation of scale, changing from the order of cm 3 to km 3 .The inversion scheme is typically applied to density and susceptibility models with a grid spacing of 3 x 3 x 1 km 3 (or even 5 x 5 x 2 km 3 in some large-scale models).This will yield geophysical properties averaged over a volume of 9 to 50 km 3 .In this case, the rock-physics inversion cannot and should not capture the wide spread of the observed values in the rock samples, neither in terms of the crustal properties (density, susceptibility) nor in inverted RHP values.Therefore, the rock-physics models have been designed a priori to avoid very high RHP values (above 5 μW/m 3 ), which are unlikely to occur in grid cells of several km 3 in size.Last but not least, a fair local geological knowledge will always be needed for interpreting the results and identifying possible false RHP predictions.

Heat production and heat flow in the greater Barents Sea
We demonstrate the full workflow on the greater Barents Sea, including LitMod3D modelling, gravity and magnetic inversion, and rock-physics inversion.
Our study area extends from the Loppa High to Novaya Zemlya, and from Finnmark to Svalbard.Input data for the study were obtained from various proprietary and public sources.The geological evolution of the area is extensively discussed in various studies (e.g., Faleide et al., 2008;Marello et al., 2013, and references therein).
The free-air gravity and magnetic data are from the GETECH compilation and NGU, respectively (Fig. 6).Regional maps of top of the crust and Moho were obtained from the work of Marello et al. (2013).The initial LAB was computed by GFZ Potsdam, and published by Klitzke et al. (2015).Seabed and regional sediment horizons were provided by Statoil Harstad.
A regional seismic velocity model was built by Statoil's depth conversion team.The horizontal sampling of all grids and models was 3 km 3 km.
First, we estimated the mantle heat flow using the LitMod3D.Parameters used in the LitMod model are shown in Table 1.We adjusted the depth to the LAB until we got a reasonable fit between modelled and observed gravity (both free-air and Bouguer), filtered geoid anomaly and bathymetry, and with the other  Koren, 2006) to the seismic velocity cube by log-linear regression, and then transferred this exponential trend to density-depth functions for the sediments (Fig. 8A).
The densities calculated from seismic velocities were then used to remove (strip off) the gravity effect of the sediments (Fig. 8B).
For the magnetic data, we assumed that the magnetic anomaly was entirely associated with the crust.The only correction needed is the reduction to pole.The inversion of the crustal gravity and magnetic anomalies was constrained by crustal geometry, using the inversion scheme described in the Appendix.The constraining data available.The depth of the final LAB varies from approximately 150 km to 190-200 km (Fig. 7A).The Moho depth varies from 29 to about 39 km in the offshore areas.Moho deeper than 40 km was modelled only onshore Norway.This causes a lateral variation in the mantle heat flow in the order of ~4 mW/ m 2 , decreasing from 21 mW/m 2 in the west to 17 mW/m 2 in the north and in the east (Fig. 7B).
Second, we performed 3D inversion of gravity and magnetics.We need to precondition the gravity data to isolate the crustal anomaly, removing the effects of the water layer, the sediments and Moho undulations (Fig. 8).
Both density and seismic velocity in sedimentary rocks are to a first order controlled by porosity.In addition, seismic velocity and density are similarly controlled by the same compaction trend.We used a 3D seismic   greater Barents Sea (Fig. 10) is about 1.5 µW/m 3 , slightly less than what we found from the onshore rock samples (Fig. 5).
To obtain the crustal contribution to the heat flow, we integrated the radiogenic heat production from the top of the crystalline crust to the top of the lower crust (Fig. 8A).By adding the mantle contribution (Fig. 7B) and the heat production in the sediments, we obtain an estimate of the surface heat flow (Fig. 11).The basement highs stand out with elevated heat production and heat flow, possibly because of thicker upper crust and more density anomalies from the inversion were then added to the background density model (i.e., the model used in the gravity stripping) to obtain absolute densities.For magnetic susceptibility, the background model is zero.
Third, we computed radiogenic heat production in the crust by Bayesian rock-physics inversion of the density and susceptibility models (Fig. 9).Radiogenic elements are not evenly distributed in the crust, but tend to be concentrated in the upper part.Gravity and magnetic data are, however, not able to resolve this variation in detail.The computed average heat production in the 11).The heat-flow analysis in the wells is based on Bullard's plots (Beardsmore & Cull, 2001), and is using the relationship between thermal conductivity and seismic velocity presented by Duffaut et al. (2017).An example from the Johan Castberg area is shown in Fig. 12.
Given the present-day surface heat flow, we can compute the steady-state approximation to the subsurface felsic compositions.In this regional study, we have not removed the effect of salt diapirs from the gravity data.Therefore, the heat production and heat flow are erroneous in the Nordkapp Basin and Tiddly Bank Basin.
We can compare and, if necessary, calibrate and adjust the surface heat-flow maps to heat flow computed from borehole temperatures (indicated by the circles in Fig. Figure 12.Surface heat flow (q 0 ), seismic velocity (v P ), thermal conductivity (k), temperature gradient (T g ), and temperature (T), computed from borehole temperatures (green dots), acoustic log and gamma-ray log, using the rock-physics relation of Duffaut et al., (2017).The example is from the Johan Castberg area.areal coverage.Seismic velocity models of the crust can also be utilised if available.
The proposed workflow gives a present-day snapshot of the heat flow.It is applicable to present-day temperature and maximum paleo-temperature prediction in the steady-state approximation, given that the heat flow has not changed since the maximum burial.This is an acceptable assumption for the Barents Sea.To obtain the heat-flow history for time-dependent thermal modelling, techniques from kinematic restoration and tectonic modelling must be applied, e.g., using the McKenzie (1978) model of crustal stretching.In this case, the estimated heat production and present-day heat flow can be used to constrain the kinematic restoration.
We demonstrated the proposed inversion method in two different case studies from onshore Norway, using data from rock samples, and from the Barents Sea using gravity and magnetic data.For the Barents Sea crust, we found an average radiogenic heat production of approximately 1.5 µW/m 3.This is slightly less than for the onshore rock samples.Also, in the Barents Sea case, we do not find the very high values (A > 5 µW/m 3 ) that we estimated for some of the onshore samples.
The onshore data are measured on small near-surface rock samples, which are more likely to give very high or very low values.However, the estimates from the Barents Sea are averages over the entire upper and middle crust, which is the best estimate that can be resolved with gravity and magnetic data.Rock samples may encounter very high local heat production, which is not found on the kilometre scale, on which the gravity and magnetic inversions work.Our estimated heat-flow values in the Barents Sea are within the range reported by the NGU HeatBar Project (Pascal et al., 2010).
temperature.Combining Fourier's law (Equation 1) with the relationship between thermal conductivity and seismic velocity presented by Duffaut et al. (2017), we can compute the temperature of target formations from seismic time horizons (Hokstad, 2014).Moreover, using a map of net erosion (Henriksen et al., 2011), we can estimate the temperature at maximum burial.As an example, we show the computed present-day and maximum paleo-temperature at the Top Klappmyss Formation (Fig. 13).

Discussion and conclusions
We have presented and demonstrated a novel methodology for estimation of radiogenic heat production in the crust, and basal and surface heat flow.
The results give valuable input to basin modelling.
The mantle contribution to the heat flow is mainly determined by the depth to the LAB and Moho.This component is obtained using the LitMod3D modelling software, constrained by elevation, gravity data and filtered geoid anomaly.In some cases, the depth to the LAB is uncertain, and we need to work with two or more scenarios for upper and lower bounds.
The crustal contribution to the heat flow is computed by multi-geophysical inversion, including Bayesian rockphysics inversion.The inversion can be performed with only one input dataset, e.g., gravity, but the use of at least two different types of geophysical data reduces the uncertainty in estimated radiogenic heat production and heat flow.Most important for regional exploration are gravity and magnetic data, which are available with large

Appendix: Gravity and magnetic inversion scheme
After appropriate preconditioning, the vertical component of the gravity response ∆ O density anomaly ∆ is given by Newton's law of gravity where  R is the gravity constant.Boldface symbols denote vectors or matrices.The compo of the gravity tensor ∆ /X are given by the derivatives of the gravity field Here, we will only consider the vertical derivative of the vertical component From Equations A-1 and A-3, it is clear that gravity data are linear in the densities (but n the positions of density anomalies), and lead to a linear inverse problem.Hence, discretising integrals above, we can write the forward problem as where  is the (Gaussian) noise, and the components of the model vector  are  X = ∆  X .
For standard gravity data ∆ O the elements of the data vector  and matrix  are where ∆ X is a finite volume element.For tensor gravity ∆ OO , the elements of the data vect and matrix  are

𝐴𝐴
From Equations A-1 and A-3, it is clear that gravity data are linear in the densities (but n the positions of density anomalies), and lead to a linear inverse problem.Hence, discretising integrals above, we can write the forward problem as where  is the (Gaussian) noise, and the components of the model vector  are  X = ∆  X .
For standard gravity data ∆ O the elements of the data vector  and matrix  are where ∆ X is a finite volume element.For tensor gravity ∆ OO , the elements of the data vect and matrix  are

𝐴𝐴
From Equations A-1 and A-3, it is clear that gravity data are linear in the densities (but n the positions of density anomalies), and lead to a linear inverse problem.Hence, discretising integrals above, we can write the forward problem as where  is the (Gaussian) noise, and the components of the model vector  are  X = ∆  X .
For standard gravity data ∆ O the elements of the data vector  and matrix  are where ∆ X is a finite volume element.For tensor gravity ∆ OO , the elements of the data vect and matrix  are

Appendix: Gravity and magnetic inversion scheme
After appropriate preconditioning, the vertical component of the gravity response ∆ O o density anomaly ∆ is given by Newton's law of gravity where  R is the gravity constant.Boldface symbols denote vectors or matrices.The compone of the gravity tensor ∆ /X are given by the derivatives of the gravity field Here, we will only consider the vertical derivative of the vertical component Appendix: Gravity and magnetic inversion scheme After appropriate preconditioning, the vertical component of the gravity response ∆ O density anomaly ∆ is given by Newton's law of gravity where  R is the gravity constant.Boldface symbols denote vectors or matrices.The compo of the gravity tensor ∆ /X are given by the derivatives of the gravity field Here, we will only consider the vertical derivative of the vertical component Appendix: Gravity and magnetic inversion scheme After appropriate preconditioning, the vertical component of the gravity response ∆ O density anomaly ∆ is given by Newton's law of gravity where  R is the gravity constant.Boldface symbols denote vectors or matrices.The comp of the gravity tensor ∆ /X are given by the derivatives of the gravity field Here, we will only consider the vertical derivative of the vertical component all number that makes 1/ X finite everywhere (typically  = 0.001).In the ion A-10 we then choose diagonal covariance matrices as Assuming reduction to pole, and neglecting remanent magnetisation, the vertical componen the magnetic data is given by where  =  is the magnitude of the magnetic moment The magnetic field is given by Introducing a susceptibility anomaly ∆, Equation A-19 can be rewritten as where  is the background reference magnetic field (typically  ≈ 50000 nT).
Magnetic data are linear in the susceptibilities, with the same Green's function as the gra gradient  OO .This leads to the same linear inverse problem as for gravity discussed ab Discretising the integral above, we obtain a linear system of equations with elements ed horizons to constrain geometry, and inverting only for properties.Given top ons  h ,  and  n ,  , we define weight functions all number that makes 1/ X finite everywhere (typically  = 0.001).In the ion A-10 we then choose diagonal covariance matrices as all number that makes 1/ X finite everywhere (typically  = 0.001).In the ion A-10 we then choose diagonal covariance matrices as all number that makes 1/ X finite everywhere (typically  = 0.001).In the ion A-10 we then choose diagonal covariance matrices as

25
Assuming reduction to pole, and neglecting remanent magnetisation, the vertical compo the magnetic data is given by where  =  is the magnitude of the magnetic moment
Introducing a susceptibility anomaly ∆, Equation A-19 can be rewritten as where  is the background reference magnetic field (typically  ≈ 50000 nT).
Magnetic data are linear in the susceptibilities, with the same Green's function as the gradient  OO .This leads to the same linear inverse problem as for gravity discussed Discretising the integral above, we obtain a linear system of equations with elements where  is the background reference magnetic field (typically  ≈ 50000 nT).
Magnetic data are linear in the susceptibilities, with the same Green's function as the gradient  OO .This leads to the same linear inverse problem as for gravity discussed Discretising the integral above, we obtain a linear system of equations with elements where  is the background reference magnetic field (typically  ≈ 50000 nT).
Magnetic data are linear in the susceptibilities, with the same Green's function as the gra gradient  OO .This leads to the same linear inverse problem as for gravity discussed ab Discretising the integral above, we obtain a linear system of equations with elements where  is the background reference magnetic field (typically  ≈ 50000 nT).
Magnetic data are linear in the susceptibilities, with the same Green's function as the g gradient  OO .This leads to the same linear inverse problem as for gravity discussed Discretising the integral above, we obtain a linear system of equations with elements

Figure 1 .
Figure 1.Bayesian network for estimation of radiogenic heat production A from gravity g z , magnetic T z , and (optionally) seismic data P.The red frame indicates the part computed by Bayesian inversion of susceptibility χ, density ρ, and (optionally) P-wave velocity v P .
is the normalisation factor, and   is the prior distribution for the radiogenic production.Different lithology-dependent priors can be used if we have relevant knowled the nature of the crust (oceanic or continental).For the problem considered here,  / = { where  is density and  is magnetic susceptibility.Assuming Gaussian errors, each o likelihood functions   / | in the product on the right-hand side of Equation (2), c written as

Figure 3 .
Figure 3. (A) Observed RHP in µW/m 3 computed from measured uranium, thorium and potassium contents.(B) RHP in µW/m 3 computed by Bayesian rock-physics inversion of measured density and susceptibility.The RHP is plotted on the geographical origin of each sample.

Figure 4 .
Figure 4. Cross-plot of observed vs. inverted RHP averaged per lithology type for the samples shown in Fig. 3.

Figure 5 .
Figure 5. (A) Cross-plot of measured density vs. magnetic susceptibility from approximately 18,000 rock samples acquired onshore Norway.The colours indicate crustal lithology.(B) RHP computed by Bayesian rockphysics inversion of the density and susceptibility data in A, colour-coded by RHP, and plotted on the geographical origin of each sample.Data from NGU(Olesen et al., 2010).

Figure 6 .
Figure 6.(A) Free-air gravity.(B) Total magnetic intensity, reduced to pole.The blue line is the border between Norway and Russia.

Figure 8 .
Figure 8. (A) Vertical section from the background density model.(B) Crustal gravity anomaly after stripping off the effects of the water layer, sediments and Moho undulations.The red line shows the location of the vertical section in A

Figure 9 .
Figure 9. Vertical section of input data and output models from the multi-geophysical inversion.(A) Mantle heat flow and basal heat flow in mW/m 2 .(B) RHP in µW/m 3 .(C) Inverted density.(D) Inverted susceptibility.Top of crust, top of lower crust and Moho are indicated with black lines.The location of the vertical section is indicated by the red line in Figure 8B.

Figure 10 .
Figure 10.Average RHP in µW/m 3 computed by inversion in the Barents Sea crust.

Figure 11 .
Figure 11.Surface heat-flow map from multi-geophysical inversion.Circles show heat flow estimated from measured (and Horner corrected) borehole temperatures, acoustic logs and gamma-ray logs.

Figure 13 .
Figure 13.(A) Present-day temperature at Top Klappmyss Fm. (B) Max paleo-temperature at Top Klappmyss Fm.Temperatures in °C were computed using the method of Hokstad (2014).Inputs to the calculation are the surface heat-flow map in Fig 11, the interpreted seismic time horizon, and for max paleo-temperature the net erosion map presented byHenriksen et al. (2011).
T F  U  T , where  =  is the magnitude of the magnetic moment = .The magnetic field is given by =   +  =  1 +  .Introducing a susceptibility anomaly ∆, Equation A-19 can be rewritten as =  is the magnitude of the magnetic moment = . −The magnetic field is given by =   +  =  1 +  . −Introducing a susceptibility anomaly ∆, Equation A-19 can be rewritten as T F  U  T ,  where  =  is the magnitude of the magnetic moment = .The magnetic field is given by =   +  =  1 +  .Introducing a susceptibility anomaly ∆, Equation A-19 can be rewritten as

Table 1 .
Hasterok (2010) and the forward modelling with LitMod 3D.Pressure/depth-dependent densities (*) of the sediments are constrained by the seismic velocity model used for the depth-density functions.P/T/C dependent (C stands for composition) densities (**) are given as ranges.The RHP values in the sediments are constrained by in-house data; the mantle values by Hasterok 2010; the other values are constrained byHasterok (2010) and Vilà et al. (2010)and are partly a result of the modelling.
T F  U  T . finite volume element.For tensor gravity ∆ OO , the elements of the data vector  Oldenburg (1998) is commonly used.Here, we propose a different approach, ed horizons to constrain geometry, and inverting only for properties.Given top ons  h ,  and  n ,  , we define weight functions