Constraining the geotherm beneath the British Isles from Bayesian inversion of Curie depth: integrated modelling of magnetic, geothermal, and seismic data

Curie depth offers a valuable constraint on the thermal structure of the lithosphere, based on its interpretation as the depth to 580 ^∘C, but current methods underestimate the range of uncertainty. We formulate the estimation of Curie depth within a Bayesian framework to quantify its uncertainty across the British Isles. Uncertainty increases exponentially with Curie depth but this can be moderated by increasing the size of the spatial window taken from the magnetic anomaly. The choice of window size needed to resolve the magnetic thickness is often ambiguous but, based on our chosen spectral method, we determine that significant gains in precision can be obtained with window sizes 15–30 times larger than the deepest magnetic source. Our Curie depth map of the British Isles includes a combination of window sizes: smaller windows are used where the magnetic base is shallow to resolve small-scale features, and larger window sizes are used where the magnetic base is deep in order to improve precision. On average, the Curie depth increases from Laurentian crust (22.2±5.3 km) to Avalonian crust (31.2±9.2 km). The temperature distribution in the crust, and associated uncertainty, was simulated from the ensemble of Curie depth realizations assigned to a lower thermal boundary condition of a crustal model (sedimentary thickness, Moho depth, heat production, thermal conductivity), constructed from various geophysical and geochemical datasets. The uncertainty in the simulated heat flow field substantially increases from ±10 mW m⁻² for shallow Curie depths at ∼15 km to ±80 mW m⁻² for Curie depths >40 km. Surface heat flow observations are concordant with the simulated heat flow field except in regions that contain igneous bodies. Heat flow data within large batholiths in the British Isles exceed the simulated heat flow by ∼25 mW m⁻² as a result of their high rates of heat production (4–6 µW m⁻³). Conversely, heat refraction around thermally resistive mafic volcanics and thick sedimentary layers induce a negative heat flow misfit of a similar magnitude. A northward thinning of the lithosphere is supported by shallower Curie depths on the northern side of the Iapetus Suture, which separates Laurentian and Avalonian terranes. Cenozoic volcanism in Northern Britain and Ireland has previously been attributed to a lateral branch of the proto-Icelandic mantle plume. Our results show that high surface heat flow (>90 mW m⁻²) and shallow Curie depth (∼15 km) occur within the same region, which supports the hypothesis that lithospheric thinning occurred due to the influence of a mantle plume. The fact that the uncertainty is only ±3–8 km in this region demonstrates that Curie depths are more reliable in hotter regions of the crust where the magnetic base is shallow.

1 Introduction

Magnetic data are, along with gravity anomalies, the most commonly available geophysical observables for subsurface imaging of the Earth. Surveys of the magnetic anomaly capture, among other features, the contribution of various magnetic minerals to the surface field recorded by magnetometers. Ferromagnetic minerals retain their magnetism until they reach their Curie temperature. Magnetite is the most prevalent magnetic mineral, in terms of susceptibility and quantity, and has a Curie temperature of 580 ^∘C. Therefore, the Curie depth is often interpreted as the depth to the 580 ^∘C isotherm. For this reason, Curie depth offers a valuable constraint to estimate the geothermal heat flow.

Various methods have been proposed to estimate Curie depth from a spectral analysis of the magnetic anomaly, which have been successfully applied to many regions of the Earth, but often under-represent the degree of uncertainty within each Curie depth estimate. We propose a Bayesian approach where Curie depth is expressed in probabilistic terms. We apply this methodology to the British Isles to quantify the degree of uncertainty in Curie depth and use this as a boundary condition to compute the temperature distribution in the crust, integrating a stratified geological model constrained by heat production, seismic, and surface heat flow data. We validate our prediction of geothermal heat flow across the British Isles against surface heat flow data to examine the variations and controls on crustal heat flow.

1.1 Geological and tectonic context

The British Isles comprise Laurentian and Avalonian terranes that were brought together at the closure of the Iapetus Ocean at 420 Ma (Chew and Strachan, 2014) (Fig. 1). The Iapetus Suture Zone (ISZ) is an area of deformation along which continental accretion took place during the Caledonian orogeny (Soper et al., 1992). Seismic studies identify the ISZ as a thickened mid-crust along a 60 km wide zone beneath central Ireland and northern England (Klemperer et al., 1991; Brewer et al., 1983). The major exposed Caledonian granite batholiths were emplaced during widespread magmatism from the end of the Silurian to the Devonian (Brown et al., 2008). Extension associated with sea-floor spreading in the North Atlantic during the Cretaceous formed sedimentary basins off the west coast of Ireland (e.g. O'Reilly and Griffin, 2010). Widespread extrusion of flood basalts during the Cenozoic, possibly associated with a branch of the Icelandic mantle plume (Landes et al., 2007), formed the North Atlantic Igneous Province between Northern Ireland and Scotland.

The V_p structure of the British Isles generally supports a uniformly thick crust in Avalonia of 32–35 km (Landes et al., 2005; Davis et al., 2012), which thins north of the ISZ from 35 to 25 km in northwest Scotland and Ireland (Tomlinson et al., 2006; Licciardi et al., 2014). The Laurentian side of the British Isles has experienced significant igneous activity, initially during Caledonian emplacement of granite batholiths and later due to extrusion of mafic volcanics. High radiogenic heat production rates (4 to 8 µW m⁻³) have been inferred in the Laurentian terrane of Ireland and Scotland (Webb et al., 1987; Mather et al., 2018). In contrast, the comparatively fewer granitic bodies outcropping in the Avalonian crust are associated with lower heat production rates (3.5–5.5 µW m⁻³; Webb et al., 1987).

Figure 1 Igneous outcrop and major geological structures across the British Isles. Surface heat flow data are compiled from Mather et al. (2018) and Pollack et al. (1993).

2 Methodology

We apply probabilistic techniques to quantify the uncertainty in Curie depth (interpreted as the 580 ^∘C isotherm) estimated from magnetic anomaly and the associated surface heat flow. In particular, we cast the computation of Curie depth using spectral methods within a Bayesian framework. The Curie depth is then used as the lower boundary condition of a 3-D Cartesian mesh, over which we solve the steady-state heat equation, for a given set of thermal parameters, in order to estimate the surface heat flow and its uncertainty.

2.1 Spectral methods to compute Curie depth

Most methods to estimate Curie depth relate the spectrum of magnetic anomalies and the depth of magnetic sources by transforming the data into the Fourier domain, which is computed over a square window of the magnetic anomaly. Depth to the bottom of magnetic sources is estimated from the slope of the radial power spectrum, Φ(k),

$$\begin{array}{}\text{(1)}& \mathrm{\Phi }\left(k\right)\propto k^{-\mathit{\beta }},\end{array}$$

where k is the wave number in Cartesian space, $k=\left|k\right|$ is its norm, and β is the fractal parameter of magnetization. Originally, Tanaka et al. (1999) analysed the slope of the radial power spectrum at high and low wave numbers to estimate the top and bottom of magnetic sources, respectively. Their method has been applied to many regions worldwide with small modifications (Li et al., 2009; Li, 2011; Salem et al., 2014; Salazar et al., 2017; Martos et al., 2017). In this method of Curie depth estimation, the wave number segments over which to calculate the depth of magnetic sources vary in different studies, but are usually <0.05 km⁻¹ (Wang and Li, 2018). However, trimming Φ to a specific range makes uncertainties difficult to quantify.

Alternatively, Bouligand et al. (2009) offered an analytic expression to the radial power spectrum that expanded on the formulation of Maus et al. (1997),

$$\begin{array}{ll}{\displaystyle }\mathrm{\Phi }\left(k\right)& {\displaystyle }=C-\mathrm{2}k{z}_{\mathrm{t}}-k\mathrm{\Delta }z-(\mathit{\beta }-\mathrm{1})\ln \left(k\right)\\ {\displaystyle }& {\displaystyle }+\left[-k\mathrm{\Delta }z+\ln \left({\displaystyle \frac{\sqrt{\mathit{\pi }}}{\mathrm{\Gamma }(\mathrm{1}+\frac{\mathit{\beta }}{\mathrm{2}})}}\left({\displaystyle \frac{\cos h\left(k\mathrm{\Delta }z\right)}{\mathrm{2}}}\mathrm{\Gamma }\left({\displaystyle \frac{\mathrm{1}+\mathit{\beta }}{\mathrm{2}}}\right)\right.\right.\right.\\ \text{(2)}& {\displaystyle }& {\displaystyle }\left.\left.\left.-K_{\frac{\mathrm{1}+\mathit{\beta }}{\mathrm{2}}}\left(k\mathrm{\Delta }z\right){\left({\displaystyle \frac{k\mathrm{\Delta }z}{\mathrm{2}}}\right)}^{\frac{\mathrm{1}+\mathit{\beta }}{\mathrm{2}}}\right)\right)\right],\end{array}$$

where the shape of the power spectrum is controlled by four variables: β – a fractal parameter, z_t – the top of magnetic sources, Δz – the thickness of the magnetic layer, and C – a field constant. The Curie depth is found by summation of z_t and Δz. Some combination of these parameters should produce a curve that fits the radial power spectrum computed from a fast Fourier transform (FFT) of the magnetic anomaly. Li et al. (2013) demonstrated that the radial power spectrum of Maus et al. (1997) and Bouligand et al. (2009) are nearly identical, except for a constant vertical shift. Because the thickness of a buried magnetic layer depends only on the slope of the power spectrum, they yield equivalent results. Computing Curie depth in this manner offers a significant advantage over other methods because the fit between analytic and computed power spectrum can be quantified.

Figure 2 The relationship between window size and the uncertainty in Curie depth estimates (z_t+Δz) for a synthetic model. (a) Synthetic magnetic anomaly generated using the following parameters: β=3, z_t=0.305, and Δz=10; (b) radial power spectrum for different window sizes – inverted parameters for each window size are 300 km: β=2.7, z_t=0.3, Δz=8.5, $C=-\mathrm{17.6}$, 200 km: β=2.7, z_t=0.3, Δz=8.5, $C=-\mathrm{18.5}$, 100 km: β=2.8, z_t=0.3, Δz=4.8, $C=-\mathrm{19.8}$, and 50 km: β=2.9, z_t=0.3, Δz=4.9, $C=-\mathrm{21.4}$; (c) uncertainty in Curie depth estimates as a function of window size.

2.2 The inverse problem of Curie depth

Here we consider the analytic expression from Eq. (2) to compute the power spectra from the magnetic anomaly and cast it within a Bayesian framework, where information on input parameters are represented in probabilistic terms. The inverse solution is given by the a posteriori probability function, P(m|d), which is proportional to the product of the likelihood function P(d|m) and the a priori probability P(m),

$$\begin{array}{}\text{(3)}& P\left(m\right|d)\propto P(d\left|m\right)\cdot P\left(m\right).\end{array}$$

The likelihood function is the probability of reproducing the data d given a particular model m, and the a priori probability is prior knowledge about the model before assimilating the data. In this case, the model parameters correspond to the four unknowns in Eq. (2) and the data are the radial power spectrum computed from the magnetic anomaly, Φ_d; thus $P\left(m\right|d)=P(\mathit{\beta },z_{\mathrm{t}},\mathrm{\Delta }z,C\left|{\mathrm{\Phi }}_{\mathrm{d}}\right)$. The posterior probability can be evaluated through an objective function, S(m), which compares the misfit between data and prior information, m_p,

$$\begin{array}{}\text{(4)}& P\left(m\right|d)=A\exp (-S\left(m\right)),\end{array}$$

where A is a constant. The maximum a posteriori (MAP) estimate is obtained by minimizing the ℓ_p-norm objective function if data and prior information are both uncorrelated,

$$\begin{array}{}\text{(5)}& S\left(m\right)={\displaystyle \frac{\mathrm{1}}{s}}\sum _i{\displaystyle \frac{\left|g^i\right(m)-d^i{|}^s}{({\mathit{\sigma }}_{\mathrm{d}}^i{)}^s}}+{\displaystyle \frac{\mathrm{1}}{r}}\sum _j{\displaystyle \frac{\mathrm{|}{m}^j-m_{\mathrm{p}}^j{\mathrm{|}}^r}{({\mathit{\sigma }}_{\mathrm{p}}^j{)}^r}},\end{array}$$

where g is the forward operator, which is the prediction of observations from the model, m. This is simply the calculation of Φ from β, z_t, Δz, and C in Eq. (2). Gaussian probabilities are assumed for Φ (s=2, i.e. the ℓ₂-norm) and uniform priors over a large range (thus the second term in Eq. 5 disappears). The Metropolis–Hastings algorithm is implemented to sample the posterior using a “random walk” (the algorithm is detailed in Appendix A).

2.3 Resolution of model parameters: synthetic tests

The Bayesian framework we have described produces an ensemble of models that sample the posterior density function. Here we design a synthetic sensitivity test to explore the inversion parameters. As the inversion operates in the spectral domain a crucial parameter is the size of the spatial window used to transform the spatially distributed magnetic data to the Fourier domain. A synthetic magnetic anomaly was generated from random fractal noise in three dimensions with known model parameters (β=3, z_t=0.305, and Δz=10) from which spatial windows were extracted of different sizes to estimate the posterior uncertainty (Fig. 2a). Magnetic sources that reside deeper in the crust affect longer wavelengths of the magnetic field. Therefore, in order to retrieve Curie depth, the window size must be large enough to resolve low wave numbers in the radial power spectrum (Fig. 2b). Intuitively, fewer points in the low-frequency part of the radial power spectrum result in higher uncertainty; however, the relationship between window size and Curie depth uncertainty is non-linear (Fig. 2c). There are significant gains in precision for window sizes up to 200 km, and smaller improvement thereon. This non-linear relationship can be exploited to minimize the uncertainty and the computational time. Bouligand et al. (2009) estimate that window size should be 10–15 times larger than the depth extent of magnetic sources. Our results suggest that for magnetic sources in the depth range 10–40 km and uncertainties of 5 km the window size should be closer to 15–30 times the depth extent of magnetic sources.

In addition to window size, retrieval of accurate Curie depth values is also complicated by a strong correlation between β and Δz, which both control the slope of the radial power spectrum at low wave numbers (Fig. 3b). The variation in these parameters may be reduced with a priori information related to β or z_t; however, in most circumstances these parameters are completely unknown. Most studies fix β to a constant value across the entire study area; but implicit in this decision is to assume that the magnetic composition of the rocks is consistent over long wavelengths. Through casting β as an inversion variable, we retrieve a comparatively higher degree of uncertainty because our method propagates all of the errors associated with each parameter within a Bayesian framework. If this parameter is not taken into account then it is likely that all other parameters (z_t, Δz, and C) will be biased.

Figure 3 Ensemble of radial power spectra from a MCMC (Markov chain Monte Carlo) simulation. (a) Radial power spectrum computed over a 300 km window of a synthetic magnetic anomaly – shaded region indicates the 99th percentile of all models in the ensemble and the red curve is the MAP estimate; (b) joint posterior probability density for β and Δz – the red dot indicates MAP estimate.

2.4 Modelling geothermal heat flow

Geothermal heat flow can be determined by modelling a geotherm from the surface of the crust to the Curie depth. A geotherm is modelled by solving the steady-state heat equation,

$$\begin{array}{}\text{(6)}& \mathrm{\nabla }\left(\mathit{\lambda }\right(\mathrm{\nabla }T\left)\right)=-H,\end{array}$$

where λ is the thermal conductivity and H is the rate of heat production. Equation (6) is solved in 3-D numerically with top and bottom Dirichlet boundary conditions set to the mean annual surface temperature and Curie temperature of magnetite (12 and 580 ^∘C, respectively). The top of the model is the surface topography (Amante, 2009), and the base of the model is adjusted to the Curie depth from the inversion of magnetic anomaly as outlined in Sect. 2.1 and 2.2. To account for the decrease in thermal conductivity with temperature, λ is iteratively updated according to the relationship of Durham et al. (1987),

$$\begin{array}{}\text{(7)}& \mathit{\lambda }\left(T\right)=\mathrm{2.26}-{\displaystyle \frac{\mathrm{618.241}}{T}}+{\mathit{\lambda }}_{\mathrm{0}}\left({\displaystyle \frac{\mathrm{355.576}}{T}}-\mathrm{0.30247}\right),\end{array}$$

where λ₀ is the thermal conductivity at the surface. Surface heat flow, q_s, is calculated from the product of thermal conductivity and temperature gradient in the topmost part of the model once the temperature solution has converged after several iterations ($|T_n-T_{n-\mathrm{1}}|<{\mathrm{10}}^{-\mathrm{12}}$).

3 Results

Maps of Curie depth in the study region were taken from the magnetic anomaly of the British Isles that was extracted from the EMAG2 global compilation (Fig. 4). While its resolution of 2 arcmin (approximately 4 km) is low compared to some regional surveys, the global extent permits Curie depth estimates over very large window sizes. It is important to note that the effective resolution of this global compilation is inherited from multiple regional grids and satellite data that are spliced together to form the EMAG2 dataset. Geographical coordinates were transformed to a local Cartesian projection prior to applying the Fourier transform for the radial power spectrum to be in units of radians per kilometre. We repeated the Markov chain Monte Carlo (MCMC) simulation for 200 to 800 km windows of the magnetic anomaly at 25 km increments. A lateral rolling average of 20 km×20 km was chosen for the centroid spacing of each overlapping window. The posterior distribution is best fit by a skewed normal probability distribution function where the MAP estimate is the mode of the distribution and the uncertainty is its standard deviation.

The implemented map-generation algorithm starts at the MAP estimate for the largest window size (800 km) at every point in the model domain and iteratively reduces at 25 km increments until the standard deviation of the posterior exceeds that of the maximum window size. All other Curie depth estimates using smaller or larger window sizes for a given point are rejected. As a consequence, the inverted Curie depth map contains a combination of window sizes: smaller windows are used where the magnetic base is shallow to resolve small-scale features and larger windows are used where the magnetic base is deep in order to improve precision (Fig. 5a). While Curie depth was estimated to be 200 km, no regions were used from this window size due to the high uncertainty.

3.1 Curie depth of the British Isles and its uncertainty

A broad linear increase in the MAP estimate of Curie depth is observed from north to south of the British Isles (Fig. 6a). On average, Curie depth increases from 22.2±5.3 km in Laurentian crust to 31.2±9.2 km in Avalonian crust, concordant with a thinning of the crust north of the ISZ (Tomlinson et al., 2006; Licciardi et al., 2014). In Ireland, the Curie depth is deeper in the southwest, which points to a cooler geotherm in line with the comparatively higher electrical resistivity found in a magnetotelluric survey southwest of the ISZ (Rao et al., 2014). In contrast with low electrical resistivity (for which there are various possible causes), high electrical resistivity can only be linked to dry and cold rocks (e.g. Pommier and Garnero, 2014; Khan, 2016; Fullea, 2017). A similar SW–NE pattern for shear wave velocity anomaly, with positive anomalies in the southwest, is present in a European seismic tomography model based on surface and S waveforms (Legendre et al., 2012). A European-scale adjoint waveform tomography model shows fast upper mantle velocity in the eastern English and southwestern Irish margins, whereas slow upper mantle is imaged in the S–W English peninsula (Zhu et al., 2015). Mantle seismic velocities are mainly controlled by thermal variations, with high velocities associated with low temperature and vice versa (Goes et al., 2000; Trampert, 2004; Cammarano and Romanowicz, 2007; Afonso et al., 2013; Fullea et al., 2012). Integrated geophysical–petrological modelling of gravity, elevation, and seismic data also indicate a thermal lithospheric thinning (20–30 km), which would suggest shallow Curie depth, in the north of Ireland and western Scotland (Fullea et al., 2014; Baykiev et al., 2018).

Curie depth is <20 km in Northern Ireland and Scotland, which is strongly associated with the high spatial density of plutons and mafic volcanics (Fig. 1). Curie depth at the northern edge of the Midland Platform in England is also shallow (∼20 km) despite low surface heat flow in this region. A deepening is observed south of the Variscan front in Ireland (35–40 km) but does not continue across to southern England. A distinctive increase in β occurs within England (Fig. 5b), which reflects a change in geology. Bouligand et al. (2009) noted that the fractal parameter is related to geology and may vary geographically depending on rock types of geologic structure. The fact that we observe a regional variation in β across England suggests heterogeneity at the regional scale, which can only be observed because we cast β as an inversion variable. We established that β and Δz are strongly negatively correlated (Sect. 2.3); thus, were we to fix β across the British Isles, our Curie depth estimates would be biased.

For the most part, the Curie depth we determine is concordant with previous studies. The global reference model of Li et al. (2017) is also based on the same EMAG2 magnetic anomaly and follows the similar pattern of variation we observe in the British Isles, albeit at lower resolution. Their method is adapted from Tanaka et al. (1999) and does not consider the estimated uncertainty of Curie depth. In Baykiev et al. (2018), the Curie depth (without uncertainty) is determined from forward modelling gravity, seismic and surface elevation data, and the vertically integrated magnetic susceptibility of a simplified global model. Variations in the Curie depth are broadly similar to our results, except the amplitude of variation is substantially reduced to within 25–35 km depth.

We find that quantifying Curie depth within a Bayesian framework adds significant insight into the thermal structure of the crust. The ensemble of model simulations reveals that uncertainty correlates with Curie depth (Fig. 6b). For Curie depths of ∼15 km this equates to ±4 km uncertainty, but deeper Curie depths have a significantly higher uncertainty. This is controlled by window size: larger windows of the magnetic anomaly are required to resolve deeper magnetic thickness (Fig. 7a). On average, the Curie depth of the British Isles is between 10 and 45 km and requires 200–600 km windows of the magnetic anomaly; however, this extends to 70 km depth in the southwest of Ireland and east of England and requires up to 800 km window sizes. Uncertainty increases exponentially with Curie depth (Fig. 7b), which is due to fewer points in the low wave number portion of the radial power spectrum. The resulting effect on the posterior is an increasingly positive skew with Curie depth (Fig. 7c), i.e. a longer tail in the upper quartile range of the probability distribution. Very little skew is observed between 10 and 20 km Curie depth for window sizes ≥500 km, which indicates the posterior is a Gaussian normal distribution. The degree of skew increases at a much higher rate for smaller window sizes at Curie depths >40 km. Only the 800 km window size exhibits minor skew for magnetic sources deeper than 40 km. Higher uncertainty in the southern North Sea is compounded by the lack of magnetic data and thus very large window sizes in order to capture any sensitivity to the magnetic thickness. These results illustrate how the skew of the posterior distribution can diagnose window sizes that are insufficient to reliably estimate Curie depth.

Figure 4 Magnetic anomaly map of the British Isles, extracted from the EMAG2 version 3 global compilation (Meyer et al., 2017) and overlain with structural information adapted from Tomlinson et al. (2006). The Iapetus represents the boundary between Laurentia and Avalonia.

Figure 5 The ensemble of Curie depth estimates was computed across centroids spaced 20 km×20 km apart, as indicated by black dots in each panel. (a) Combination of window sizes used to generate the map of Curie depth (Fig. 6); each pixel illustrates the spatial resolution of our results; (b) maximum a posteriori distribution of β inverted from the magnetic anomaly and interpolated to the map resolution of 0.01^∘.

Figure 6 MAP estimate of Curie depth in the British Isles and its uncertainty; (a) Curie depth map is computed from 3900 centroids for a lateral resolution of 20 km; (b) Uncertainty is given as 1 standard deviation from the mean.

Figure 7 Analysis of Curie depth estimates. (a) The maximum Curie depth that can be resolved for a given window size; (b) the relationship of uncertainty, measured by 1 standard deviation, with respect to Curie depth; (c) the degree of skew exhibited by the posterior probability distribution as a function of Curie depth for different window sizes – skew is defined here as the mean Curie depth less the median.

3.2 Heat flow and its uncertainty

The heat flow distribution over the British Isles is determined by modelling a geotherm from the surface of the crust to the Curie depth. The top of the model is the surface topography from ETOPO2 (Amante, 2009), set to the mean annual surface temperature, and the base of the model is adjusted to the Curie depth that was determined for the British Isles. The volume between these boundaries is divided based on the layers compiled in Baykiev et al. (2018), which contains sedimentary and Moho (crust–mantle boundary) depth constrained by various geophysical sources (active seismic, receiver functions, gravity data) (Fig. 8). The crust is divided into an upper and lower layer, motivated from the well-established theory on the partitioning of heat-producing elements (e.g. Michaut et al., 2009), where an enriched upper-layer H_s (of thickness h_s) overlies a relatively depleted lower-layer H_r (of thickness h_r). The total heat production of the crust is the sum of these two layers,

$$\begin{array}{}\text{(8)}& \underset{\mathrm{crust}}{\int }{H}_{\mathrm{c}}\mathrm{d}h=H_{\mathrm{c}}{h}_{\mathrm{c}}=H_{\mathrm{s}}{h}_{\mathrm{s}}+H_{\mathrm{r}}{h}_{\mathrm{r}},\end{array}$$

where h_c is crustal thickness. The differentiation index describes the degree of partitioning in the crust,

$$\begin{array}{}\text{(9)}& D_{\mathrm{I}}={\displaystyle \frac{H_{\mathrm{s}}}{H_{\mathrm{c}}}},\end{array}$$

where higher D_I indicates that a heat flow province has undergone more significant crustal reworking of heat-producing elements to the uppermost crust. The vertical differentiation of heat-producing elements in the crust directly controls the curvature of the geotherm and, thus, depth to the 580 ^∘C isotherm. Thus, depth to the 580 ^∘C isotherm will increase if H_s≫H_r and h_s≪h_r. These values have been previously determined for the British Isles from Mather et al. (2018), where the rates of heat production are defined separately for the crust in the Laurentian and Avalonian terranes (Table 1). We split the total thickness of the crust (h_c) from Baykiev et al. (2018) into an upper (enriched) and lower (depleted) layer of variable thickness by rearranging Eq. (8) to respect the differentiation index across Laurentian and Avalonian crust (D_I=2.8 and 1.5, respectively). Therefore, the thickness of the upper (enriched) layer is

$$\begin{array}{}\text{(10)}& h_{\mathrm{s}}=h_{\mathrm{c}}{\displaystyle \frac{H_{\mathrm{c}}-H_{\mathrm{r}}}{H_{\mathrm{s}}-H_{\mathrm{r}}}},\end{array}$$

where H_s, H_r, and H_c are the rates of heat production for the upper, lower layer, and entire crust, respectively. The thermal properties assigned to each layer in our model are summarized in Table 1.

Our geothermal heat flow map offers very high resolution (approximately 4 km) compared to the spatial density of borehole measurements, and provides insights into the variation in thermal regimes across the British Isles (Fig. 10a). The optimal surface heat flow, q_s, was simulated by fitting the geometry of the lower flux boundary condition to the MAP estimate of Curie depth. Its uncertainty is determined from randomly perturbing the geometry of the lower boundary condition within the Curie depth uncertainty to build an ensemble of q_s models and taking the standard deviation (Fig. 10b). Here, we present a probabilistic estimate of surface heat flow based primarily on magnetic data and supported by geothermal observations and crustal architecture. Generally, heat flow is higher in Scotland and the northernmost parts of Ireland that are of Laurentian origin. The heat flow contrast on either side of the Iapetus Suture is linked to the crustal partitioning of heat-producing elements, but there is significant variation in q_s despite the relatively uniform Curie depth in the northern half of the British Isles. Much of this variation is due to the inclusion of sediments in the model which are thermally insulating (λ₀=2.2 W m⁻¹ K⁻¹).

Table 1 Thermal properties assigned to each layer in the British Isles, summarized from Mather et al. (2018). Units of λ₀ are W m⁻¹ K⁻¹ and H are µW m⁻³.

Figure 8 Maps of layer geometries extracted from the model of Baykiev et al. (2018), overlain with major tectonic structures. (a) Sedimentary thickness and (b) crustal thickness variations across the British Isles. The Iapetus Suture marks the boundary between Laurentia and Avalonia.

Figure 9 Schematic of each layer in the 3-D thermal model of the British Isles, subdivided from surfaces compiled in Baykiev et al. (2018) (Fig. 8). The model is divided by the Iapetus Suture that separates Laurentian and Avalonian crust. The thickness of the upper (enriched) and lower (depleted) crustal layers is determined by rearranging Eq. (8) for h_s and h_r, respectively, in Laurentian crust (H_c=1.4 µW m⁻³) and Avalonian crust (H_c=1.2 µW m⁻³).

Figure 10 Surface heat flow in the British Isles and its uncertainty, (a) Surface heat flow overlain with palaeoclimate-corrected heat flow observations compiled from Mather et al. (2018) in Ireland and Pollack et al. (1993) elsewhere – larger circles indicate a higher misfit compared to the simulated heat flow; (b) uncertainty of heat flow determined from an ensemble of model simulations with variable Curie depth, overlain with the locations of heat flow data in white. Circled regions include GG as Galway granite, PB as Palaeogene basalt, SWEB as southwest England batholith, NEB as northern England batholith, and EHB as eastern highlands batholith.

4 Discussion

Heat flow data are clustered mainly within coastlines and in some localized areas offshore. The simulated heat flow is mostly concordant with heat flow data except for locations above felsic bodies or mafic volcanics. These are not included within the 3-D model due to the unconstrained nature of their geometry. Granites can have very high rates of heat production which, when integrated across their thickness, contribute significantly to surface heat flow. This is most prevalent in the southwest England and northern England batholiths that contain rates of heat production between 3 and 5 µW m⁻³ (Webb et al., 1987) where a +25 mW m⁻² misfit is observed in comparison to the simulated heat flow. Positive misfit also occurs within the Galway granite in Ireland and the eastern highlands batholith of Scotland. In contrast, heat refracts around mafic volcanic rocks due to their low thermal conductivity. Heat refraction is most apparent in the northernmost part of Ireland where a ring of high heat flow points are scattered around the rim of Palaeogene basalts with low heat flow in its interior (58 mW m⁻²; Mather et al., 2018). While this effect is clearly visible in the British Isles, in locations where bedrock has not been mapped, either due to poor exposure or lack of geophysical information, using Curie depth to infer surface heat flow can be highly uncertain because short-wavelength features – that significantly influence the heat flow field – cannot be resolved. Where possible, surface heat flow derived from Curie depth should be constrained by crustal geometry and validated against heat flow data from boreholes.

The uncertainty of surface heat flow was estimated from spatially varying the depth of the lower boundary condition within the posterior distribution assembled in Sect. 3.1. The uncertainty of heat flow increases from ±10 mW m⁻² for shallow Curie depths of ∼15 km to ±80 mW m⁻² for Curie depths >40 km (Fig. 10b). Our method for determining surface heat flow from Curie depth is, therefore, most reliable in hot regions. Currently, uncertainty does not incorporate any variation in thermal properties assigned to each crustal layer or the geometry of each layer above the lower boundary condition. Uncertainty would likely increase if different thermal properties and layer geometries were considered in the ensemble of heat flow simulations, but they are not included here because already the uncertainty for Curie depths >30 km is more than the range of physically plausible heat flow. In general, Curie depth estimates are sensitive to the regional heat flow regime and cannot resolve anomalies that locally alter surface heat flow. These effects include granitic intrusions and hydrothermal advection among others. In spite of this, the misfit between simulated heat flow and data do not exceed 1 standard deviation of all thermal models in the ensemble (Fig. 10b). In locations of high misfit, assimilating alternate geophysical sources, such as seismic data (i.e. refraction, surface waves, or receiver functions), may offer better constraints on geothermal heat flow.

Thinning of the lithosphere towards the north of the ISZ has been associated with a branch of the Icelandic mantle plume (e.g. Landes et al., 2007; Al-Kindi et al., 2003; Kirstein and Timmerman, 2000). The lithosphere–asthenosphere boundary (LAB) has been shown to be 30 km shallower in the northern part of Ireland (Landes et al., 2007; Fullea et al., 2014; Baykiev et al., 2018), which is interpreted to be a result of the spreading head of the early Icelandic plume (Landes et al., 2007). A similar degree of thermal erosion, interpreted from seismic tomography data, appears to continue to the west coast of Scotland (Arrowsmith et al., 2005). This, in conjunction with Cenozoic volcanism and thin crust (20–27 km) in northern Britain and Ireland, suggests that this region is underlain by magmatic under-plating (Davis et al., 2012). The shallow Curie depths we estimate (∼15 km) and the associated simulated heat flow, in excess of 90 mW m⁻², lend further evidence to the fact that higher crustal temperatures are present in northern Britain and Ireland, possibly suggesting processes related to the Icelandic plume.

5 Conclusions

The Curie depth, often interpreted as the depth to 580 ^∘C, offers a valuable constraint to estimate geothermal heat flow at a large scale. Spectral methods estimate Curie depth from the slope of the radial power spectrum across a square window of the magnetic anomaly. While these methods are commonly used, they are subjective to the interpreter and underestimate the uncertainty in Curie depth estimates. We cast the analytical expression of the power spectrum formulated in Bouligand et al. (2009) within a Bayesian framework, where parameters relating to the thickness of magnetic sources are expressed in probabilistic terms.

The uncertainty in Curie depth increases rapidly with depth but can be tempered with larger window sizes. Windows of the magnetic anomaly need to be 15–30 times larger than the deepest possible magnetic base in the study area to resolve precise Curie depth estimates. In the British Isles, the Curie depth is broadly delineated by the SW–NE Iapetus Suture, which separates Laurentian and Avalonian continental blocks. The uncertainty exponentially increases with Curie depth, which is due to fewer points in the low wave number range of the radial power spectrum. Curie depths of >40 km correspond to uncertainties of ±20 km using a window size of 800 km. At these depths, the skew of the posterior distribution is used to indicate which window sizes are required to resolve Curie depths.

Surface heat flow is estimated by simulating the temperature distribution across a stratified model of the crust from the surface to the Curie depth. This 3-D model incorporates sedimentary thickness and vertically partitioned heat-producing elements. The simulated heat flow matches heat flow data to a reasonable degree of uncertainty, except where igneous bodies outcrop. Granitic intrusions can contribute significantly to surface heat flow due to their high rates of heat production, and heat refracts around mafic volcanics because of their low thermal conductivity. These results lend further support to lithospheric thinning in northern Britain and Ireland due to the presence of the proto-Icelandic mantle plume.

Appendix A: MCMC sampling

The commonly used Metropolis–Hastings algorithm samples the posterior distribution, P(m|d), using a random walk where $k=\mathrm{1},\mathrm{2},\mathrm{\dots },n$ for a Markov chain n samples long. The random walk is initialized at a random point in the prior distribution and progresses in accordance with the following algorithm:

1. Generate a proposal m^′ for the next sample by picking from the prior distribution P(m).

2. Calculate the acceptance ratio between each sample of the posterior $\mathit{\alpha }=P\left(m^{\prime }\right|d)/P(m_k\left|d\right)$.

3. Generate a random number, $\mathrm{0}\le u\le \mathrm{1}$, which will be used to decide if m^′ is accepted or rejected.

   1. accept if u≤α and set $m_{k+\mathrm{1}}=m^{\prime }$,

   2. reject if u>α and set $m_{k+\mathrm{1}}=m_k$ instead.

Numerous chains were initiated to ensure they converged to a similar solution.