Reference Models for Lithospheric Geoneutrino Signal

TL;DR

This study combines geophysical and geochemical data with Monte Carlo methods to estimate Earth's lithospheric geoneutrino signal, quantifying uncertainties and assessing model impacts on predictions for current and future detectors.

Contribution

It provides a comprehensive uncertainty analysis of geoneutrino signal predictions using multiple crustal models, highlighting the robustness of estimates across different geophysical assumptions.

Findings

Estimated total geoneutrino signal uncertainty is about 20%.

Different crustal models yield consistent signals within physical uncertainties.

Bulk crustal heat production is estimated at 7 ± 2 TW, with increased uncertainty.

Abstract

Debate continues on the amount and distribution of radioactive heat producing elements (i.e., U, Th, and K) in the Earth, with estimates for mantle heat production varying by an order of magnitude. Constraints on the bulk-silicate Earth's (BSE) radiogenic power also places constraints on overall BSE composition. Geoneutrino detection is a direct measure of the Earth's decay rate of Th and U. The geoneutrino signal has contributions from the local ($\sim$40$\%$) and global ($\sim$35$\%$) continental lithosphere and the underlying inaccessible mantle ($\sim$25$\%$). Geophysical models are combined with geochemical datasets to predict the geoneutrino signal at current and future geoneutrino detectors. We propagated uncertainties, both chemical and physical, through Monte Carlo methods. Estimated total signal uncertainties are on the order of $\sim$20$\%$, proportionally with geophysical and geochemical inputs contributing $\sim$30$\%$ and $\sim$70$\%$, respectively. We find that estimated signals, calculated using CRUST2.0, CRUST1.0, and LITHO1.0, are within physical uncertainty of each other, suggesting that the choice of underlying geophysical model will not change results significantly, but will shift the central value by up to $\sim$15$\%$, depending on the crustal model and detector location. Similarly, we see no significant difference between calculated layer abundances and bulk-crustal heat production when using these geophysical models. The bulk crustal heat production is calculated as $7 \pm2$~terrawatts, which includes an increase of 1~TW in uncertainty relative to previous studies. Future improvements, including uncertainty attribution and near-field modeling, are discussed.

Tables

Table 1. Table 1: Physical properties of the layers of the Bulk-Silicate Earth (BSE) modeled using the geophysical models CRUST2.0, CRUST1.0, and LITHO1.0. The density reported is the mean between the three models with variation between model layers of < < 2.5% (except 10.5% variation for lower crust). BSE mass is from a combination of PREM Dziewonski  \BBA Anderson ( \APACyear 1981 ) derived mass of the DM and EM combined with the lithosphere mass. * Thicknesses are un-weighted of the 1 × \times 1 degree data unlike the area-weighted average reported in Figure  4 .

CRUST2.0

CRUST1.0

LITHO1.0

Layer

$\sim \rho$ (g/cm3)

d (km)*

M (1021 kg)

d (km)*

M (1021 kg)

d (km)*

M (1021 kg)

CC

Sed

2.2

2.18 $\pm$2.1

0.8 $\pm$0.1

1.64 $\pm$2.2

0.7 $\pm$0.1

1.54 $\pm$2.0

0.7 $\pm$0.1

UC

2.75

11.60 $\pm$3.9

7.0 $\pm$0.9

11.71 $\pm$4.0

6.3 $\pm$0.8

12.79 $\pm$4.1

6.9 $\pm$0.8

MC

2.84

11.18 $\pm$3.4

7.2 $\pm$0.9

11.57 $\pm$3.0

6.4 $\pm$0.8

13.06 $\pm$3.8

7.3 $\pm$0.9

LC

3.02

9.93 $\pm$2.9

6.7 $\pm$0.8

10.73 $\pm$2.7

6.2 $\pm$0.8

12.22 $\pm$3.6

7.2 $\pm$0.9

Bulk CC

2.9

34.25 $\pm$8.8

21.8 $\pm$2.6

35.53 $\pm$7.6

19.6 $\pm$2.4

39.60 $\pm$9.1

22.2 $\pm$2.6

OC

Sed

1.9

1.86 $\pm$0.2

0.3 $\pm$0.0

1.90 $\pm$0.1

0.4 $\pm$0.1

1.90 $\pm$0.1

0.4 $\pm$0.1

C

2.9

6.80 $\pm$1.5

5.6 $\pm$0.7

7.82 $\pm$2.9

7.1 $\pm$0.9

10.46 $\pm$4.3

9.2 $\pm$1.1

Mantle

LM

3.34

139.9 $\pm$75

102.8 $\pm$53

139.5 $\pm$75

88 $\pm$43

114.3 $\pm$82

63 $\pm$7.5

DM

4.4

1966

3149.9

1966

3168.5

1987 $\pm$84

3187.3

EM

5.4

750

754

750

754

750

754

BSE

4.45

2891

4035

2891

4038

2891

4036

Table 2. Table 2: Calculated abundances, element ratios, and power for the upper crust, middle crust, lower crust, and bulk-continental crust. Upper crust abundances are same in all models, but power varies as mass of crust varies. Reported uncertainties are 1 σ 𝜎 \sigma . Bulk continental crust values are weighted by mass of each layer and include the sediment contribution (see Supplementary Information for sediment values). K/U and Th/U are mass ratios. a = Values from \citeA rudnick2014. b = Derived from V P subscript 𝑉 𝑃 V_{P} –(U,Th,K) relationship of amphibolites. c = Derived from V P subscript 𝑉 𝑃 V_{P} –(U,Th,K) relationship of granulites.

		CRUST2.0	CRUST1.0	LITHO1.0
Upper Crust a	U ($\mu$g/g)		2.7 $\pm$0.60
	Th ($\mu$g/g)		10.5 $\pm$1.0
	K (wt%)		2.32 $\pm$0.19
	K/U	8,900 ${}_{-1,300}{}^{}{}_{}{}^{+1,600}$	8,900 ${}_{-1,300}{}^{}{}_{}{}^{+1,600}$	8,900 ${}_{+1,300}{}^{}{}_{}{}^{+1,600}$
	Th/U	4.0 ${}_{-0.6}{}^{}{}_{}{}^{+0.7}$	4.0 ${}_{-0.6}{}^{}{}_{}{}^{+0.6}$	4.0 ${}_{-0.6}{}^{}{}_{}{}^{+0.7}$
	P (TW)	4.2 ${}_{-0.8}{}^{}{}_{}{}^{+0.9}$	3.8 ${}_{-0.7}{}^{}{}_{}{}^{+0.8}$	4.1 ${}_{-0.7}{}^{}{}_{}{}^{+0.9}$
Middle Crust b	U	0.82 ${}_{-0.42}{}^{}{}_{}{}^{+0.87}$	0.94 ${}_{-0.47}{}^{}{}_{}{}^{+0.95}$	0.84 ${}_{-0.43}{}^{}{}_{}{}^{+0.87}$
	Th	3.64 ${}_{-2.30}{}^{}{}_{}{}^{+6.27}$	4.62 ${}_{-2.80}{}^{}{}_{}{}^{+7.11}$	3.88 ${}_{-2.36}{}^{}{}_{}{}^{+6.03}$
	K	1.42 ${}_{-0.73}{}^{}{}_{}{}^{+1.52}$	1.69 ${}_{-0.84}{}^{}{}_{}{}^{+1.66}$	1.47 ${}_{-0.74}{}^{}{}_{}{}^{+1.47}$
	K/U	14,700 $\pm$2,200	15,300 $\pm$2,000	14,900 $\pm$2,000
	Th/U	4.6 $\pm$1.4	5.0 $\pm$1.2	4.7 $\pm$1.2
	P	1.4 ${}_{-0.8}{}^{}{}_{}{}^{+1.8}$	1.5 ${}_{-0.8}{}^{}{}_{}{}^{+1.7}$	1.5 ${}_{-0.8}{}^{}{}_{}{}^{+1.7}$
Lower Crust c	U	0.15 ${}_{-0.09}{}^{}{}_{}{}^{+0.22}$	0.19 ${}_{-0.11}{}^{}{}_{}{}^{+0.27}$	0.17 ${}_{-0.10}{}^{}{}_{}{}^{+0.24}$
	Th	0.81 ${}_{-0.57}{}^{}{}_{}{}^{+1.94}$	1.14 ${}_{-0.82}{}^{}{}_{}{}^{+2.89}$	0.95 ${}_{-0.67}{}^{}{}_{}{}^{+2.27}$
	K	0.70 ${}_{-0.37}{}^{}{}_{}{}^{+0.80}$	0.91 ${}_{-0.50}{}^{}{}_{}{}^{+1.12}$	0.78 ${}_{-0.42}{}^{}{}_{}{}^{+0.90}$
	K/U	39,700 $\pm$10,600	42,000 $\pm$10,400	40,700 $\pm$10,400
	Th/U	5.3${}_{-1.7}{}^{}{}_{}{}^{+2.5}$	6.0 ${}_{-2.1}{}^{}{}_{}{}^{+3.1}$	5.6 ${}_{-1.8}{}^{}{}_{}{}^{+2.6}$
	P	0.4${}_{-0.2}{}^{}{}_{}{}^{+0.6}$	0.5${}_{-0.3}{}^{}{}_{}{}^{+0.8}$	0.5${}_{-0.3}{}^{}{}_{}{}^{+0.7}$
Bulk Continental Crust	U	1.31${}_{-0.28}{}^{}{}_{}{}^{+0.35}$	1.35${}_{-0.29}{}^{}{}_{}{}^{+0.36}$	1.29${}_{-0.27}{}^{}{}_{}{}^{+0.35}$
	Th	5.81${}_{-1.47}{}^{}{}_{}{}^{+1.96}$	6.25${}_{-1.66}{}^{}{}_{}{}^{+2.26}$	5.77${}_{-1.48}{}^{}{}_{}{}^{+1.99}$
	K	1.81${}_{-0.39}{}^{}{}_{}{}^{+0.51}$	1.96${}_{-0.44}{}^{}{}_{}{}^{+0.57}$	1.81${}_{-0.40}{}^{}{}_{}{}^{+0.52}$
	K/U	11,500${}_{-1,600}{}^{}{}_{}{}^{+1,800}$	12,100${}_{-1,700}{}^{}{}_{}{}^{+2,000}$	11,800${}_{+1,700}{}^{}{}_{}{}^{+2,000}$
	Th/U	4.4${}_{-0.6}{}^{}{}_{}{}^{+0.7}$	4.6${}_{-0.7}{}^{}{}_{}{}^{+0.8}$	4.5${}_{-0.7}{}^{}{}_{}{}^{+0.8}$
	P	7.0${}_{-1.6}{}^{}{}_{}{}^{+2.1}$	6.6${}_{-1.6}{}^{}{}_{}{}^{+2.0}$	7.0${}_{-1.6}{}^{}{}_{}{}^{+2.1}$

Table 3. Table 3: Variables used in the calculation and propagation of the geoneutrino signal.

Symbol	Description	Units
$\frac{dN(E_{{\overline{\nu }}_e},\overrightarrow{r})}{d{E}_{{\overline{\nu }}_e}}$	${\overline{\nu }}_e$ detection spectrum	${\overline{\nu }}_e$
$ϵ$	${10}^{32}$ proton $\times$ 3.154$\times {10}^7$ s $\times$ 100%	proton$\times$s
$N_A$	Avagodro’s number	$\frac{atom}{mol}$
$\lambda$	Decay constant	$\frac{decay}{s\times atom}$
$\mu$	Atomic mass in kg	$\frac{kg}{mol}$
${\sigma }_P\left(E_{{\overline{\nu }}_e}\right)$	${\overline{\nu }}_e$ cross-section in $m^2$ (function of $E_{{\overline{\nu }}_e}$)	$\frac{m^2}{proton}$
$\frac{dn\left(E_{{\overline{\nu }}_e}\right)}{d{E}_{{\overline{\nu }}_e}}$	${\overline{\nu }}_e$ emission spectrum	$\frac{{\overline{\nu }}_e}{decay}$
$P_{ee}(E_{{\overline{\nu }}_e},|\overrightarrow{r}-\overrightarrow{r}\prime \left|\right)$	Oscillation probability (function of $E_{{\overline{\nu }}_e}$)	unit-less
$a(\overrightarrow{r}\prime )$	Abundance of radionuclide in cell	$\frac{kg}{kg}$
$\rho (\overrightarrow{r}\prime )$	Density of rock in cell	$\frac{kg}{m^3}$
$|\overrightarrow{r}-\overrightarrow{r}\prime |$	Distance from cell to detector	m

Table 4. Table 4: Continental crust (top row) and total (bottom row) signal (in TNU) using each detector location from each geophysical model. ∣ Δ ∣ delimited-∣∣ Δ \mid\Delta\mid is the maximum difference between the largest and smallest signal’s central value estimated at each detector, with ∣ Δ ∣ delimited-∣∣ Δ \mid\Delta\mid = ∣ ∣ \mid (x-y)/((x+y)/2) ∣ × \mid\times 100%. The overlapping coefficient refers to the degree of overlap of the two distributions listed at the head of each column (from 0 to 100% overlap; eqn. 5 ), with L1 = LITHO1.0, C1 = CRUST1.0, and C2 = CRUST2.0. Values of ∣ Δ ∣ delimited-∣∣ Δ \mid\Delta\mid ( > > 10) and overlapping coefficient ( < < 80) which show large differences to the other models are highlighted in red.

		Geoneutrino Signal (TNU)				Overlapping Coefficient (%)
Detector		CRUST2	CRUST1	LITHO1	$\mid 𝚫\mid$(%)	L1,C1	L1,C2	C1,C2
KamLAND	Bulk CC	22.7${}_{-4.7}{}^{}{}_{}{}^{+5.9}$	24.2${}_{-5.2}{}^{}{}_{}{}^{+6.7}$	26.4${}_{-5.6}{}^{}{}_{}{}^{+7.1}$	15	85	74	89
	Total	34.7${}_{-5.0}{}^{}{}_{}{}^{+5.8}$	36.6${}_{-5.5}{}^{}{}_{}{}^{+6.5}$	37.9${}_{-5.6}{}^{}{}_{}{}^{+6.6}$	9	92	78	87
Borexino		30.5${}_{-6.4}{}^{}{}_{}{}^{+8.1}$	29.9${}_{-6.3}{}^{}{}_{}{}^{+8.0}$	30.5${}_{-6.2}{}^{}{}_{}{}^{+7.7}$	2	96	98	96
		43.2${}_{-6.7}{}^{}{}_{}{}^{+8.0}$	42.9${}_{-6.6}{}^{}{}_{}{}^{+7.9}$	42.1${}_{-6.2}{}^{}{}_{}{}^{+7.3}$	1	95	94	98
SNO+		37.3${}_{-8.0}{}^{}{}_{}{}^{+10.2}$	32.9${}_{-7.4}{}^{}{}_{}{}^{+9.6}$	33.8${}_{-7.5}{}^{}{}_{}{}^{+9.6}$	13	95	84	80
		49.8${}_{-8.1}{}^{}{}_{}{}^{+9.7}$	45.7${}_{-7.7}{}^{}{}_{}{}^{+9.3}$	46.8${}_{-7.8}{}^{}{}_{}{}^{+9.4}$	8	95	86	82
JUNO		28.1${}_{-5.9}{}^{}{}_{}{}^{+7.5}$	28.0${}_{-6.1}{}^{}{}_{}{}^{+7.7}$	29.2${}_{-6.3}{}^{}{}_{}{}^{+8.0}$	4	93	93	99
		40.5${}_{-6.3}{}^{}{}_{}{}^{+7.4}$	40.7${}_{-6.4}{}^{}{}_{}{}^{+7.6}$	40.4${}_{-6.3}{}^{}{}_{}{}^{+7.4}$	1	98	99	99
Jinping		42.5${}_{-9.1}{}^{}{}_{}{}^{+11.5}$	47.2${}_{-10.0}{}^{}{}_{}{}^{+12.7}$	48.5${}_{-10.3}{}^{}{}_{}{}^{+13.1}$	13	96	78	83
		55.0${}_{-9.1}{}^{}{}_{}{}^{+10.9}$	59.9${}_{-10.1}{}^{}{}_{}{}^{+12.1}$	59.9${}_{-10.1}{}^{}{}_{}{}^{+12.2}$	9	100	81	81
Hawaii		2.3${}_{-0.5}{}^{}{}_{}{}^{+0.7}$	2.1${}_{-0.5}{}^{}{}_{}{}^{+0.6}$	2.3${}_{-0.5}{}^{}{}_{}{}^{+0.7}$	6	90	93	83
		12.9${}_{-2.3}{}^{}{}_{}{}^{+2.8}$	12.9${}_{-2.3}{}^{}{}_{}{}^{+2.8}$	12.8${}_{-2.3}{}^{}{}_{}{}^{+2.7}$	1	98	99	99

Equations

$$V_{model}=V_{f}f+V_{m}m$$

$$1=f+m$$

$$a_{output}=a_{f}f+a_{m}m$$

$$\frac{dN(E_{{\bar{\nu }}_e},\vec{r})}{d{E}_{{\bar{\nu }}_e}}=ϵ\frac{N_A\lambda }{\mu }{\sigma }_P(E_{{\bar{\nu }}_e})\frac{dn(E_{{\bar{\nu }}_e})}{d{E}_{{\bar{\nu }}_e}}{\int }_{\text{Earth}}{P}_{ee}(E_{{\bar{\nu }}_e},|\vec{r}-\vec{r}\prime |)d\vec{r}\prime \frac{a(\vec{r}\prime )\rho (\vec{r}\prime )}{4\pi |\vec{r}-\vec{r}\prime {|}^2}$$

$$OVL=\int min[f_{M_1}(S),f_{M_2}(S)]\mathrm{d}S$$

Full text

\draftfalse\journalname

arXiv

Department of Geology, University of Maryland, College Park, MD 20742, USA Department of Geophysics, Faculty of Mathematics and Physics, Charles University, Prague, Czech Republic Department of Earth Sciences and Research Center for Neutrino Science, Tohoku University, Sendai 980-8578, Japan

\correspondingauthor

S. A. Wipperfurthswipp@terpmail.umd.edu

{keypoints}

Updated geoneutrino emission model compares results using CRUST2.0, CRUST1.0, LITHO1.0 geophysical models

Observed differences between geoneutrino signals from different geophysical models is negligible

Necessary improvements to geoneutrino modeling are identified

1 Introduction

Estimated at 47 $\pm$ 2 terrawatts (TW) J\BPBIH. Davies \BBA Davies (\APACyear2010), the deep Earth’s radiant heat is primarily comprised of two sources: primordial heat remaining from planetary assembly and core formation, and radiogenic heat produced during nuclear decay of the heat-producing elements (HPE: uranium (U), thorium (Th), and potassium (K)). The radioactive isotopes of these elements — 238U (99.3% g/g of U), 235U (0.7%), 232Th (100% g/g of Th), 40K (0.012% g/g of K) — presently account for 99% of Earth’s radiogenic heat due to their long half-lives and high abundances relative to other radiogenic elements.

There are three categories of models which predict the abundance of the HPE’s in the bulk-silicate Earth (BSE; crust + mantle ($\sim$0.5% and $\sim$67% of Earth by mass, respectively)) and therefore the radiogenic heat production within the BSE. Models which predict low heat production (H) from radiogenic decay ($\approx 10$ TW) are derived from observations of isotopic similarities between Earth and enstatite chondrite Javoy (\APACyear1999); Javoy \BOthers. (\APACyear2010) or from models of early Earth collisional erosion of an HPE-enriched crust O’Neill \BBA Palme (\APACyear2008). Medium-H models ($\approx 20$ TW) are derived from combining observations from chondrites and mantle melting trends of terrestrial samples McDonough \BBA Sun (\APACyear1995); Palme \BBA O’Neill (\APACyear2014). Finally, High-H models ($\approx 30$ TW) are derived from simple parameterized mantle convection models Turcotte \BBA Schubert (\APACyear2014). There is inconclusive data to evaluate critically the veracity of each of these three BSE models, therefore there is currently not a precise understanding of the composition and thermal evolution of our planet. The Earth’s stable isotopic composition is most similar to enstatite chondrites (low-H model), and yet it falls outside of the chondritic defined endmembers in a Fe-Mg-Si plot, indicating Earth is not comprised of a single type of chondrite nor a two component mixture McDonough (\APACyear2017). Furthermore, the use of terrestrial samples and observed conservation of chondritic ratios in terrestrial samples yields a BSE composition (medium-H model) with refractory element abundances (including U and Th, but not K) double that estimated from enstatite chondrites alone. Finally, neither the low-H nor medium-H BSE models satisfy some simple parameterized convection models of the Earth which require larger amounts of radiogenic power to avoid a totally molten mantle for a significant amount of Earth history <high-H model;¿davies1980,schubert1980. Overall, the range of heat production in BSE models differ by a factor of three (10 to 30 TW). With the consideration of uncertainties and the removal of the HPE contribution from the accessible and HPE enriched continental crust <$7\pm 1$ TW;¿huang2013, these models differ by a factor of thirty in estimates of the radiogenic power in the modern mantle.

For more than a decade, particle physicists have detected and reported on the Earth’s flux of geoneutrinos — electron antineutrinos ($\bar{\nu_{e}}$) of terrestrial origin produced during $\beta^{-}$ decays <$n\rightarrow p^{+}+e^{-}+\bar{\nu_{e}}$;¿[]araki2005. The intensity of the geoneutrino flux is proportional to the concentration and sensitive to the spatial distribution of HPEs inside the Earth relative to the detector’s location. These elusive particles are exceedingly difficult to detect as they are charge-less leptons with small interaction cross sections of the weak interaction <$\sim 10^{-44}~{}cm^{2}$;¿[]vogel1999. Measurement of the geoneutrino flux requires large, underground, scintillation detectors which use the inverse beta decay (IBD) detection process ($\bar{\nu_{e}}+p^{+}\rightarrow n+e^{+}$). The IBD reaction creates two flashes of light, separated in time ($\sim$200 $\mu$s) and space ($\sim$30 cm), which uniquely classifies the event and provides an energy tag identifying the specific geoneutrino source isotope. The IBD reaction requires an $\bar{\nu_{e}}$ kinematic threshold energy of 1.806 MeV due to the larger mass of the products relative to the reactants in the production of $\bar{\nu_{e}}$’s. Consequently, only $\bar{\nu_{e}}$’s emitted by the decay of U and Th are detectable Araki \BOthers. (\APACyear2005). The detectors also register other electron antineutrinos present in the geoneutrino energy range (1.8–3.3 MeV), including reactor antineutrinos. Detectors are located at 1–2 km deep in the upper crust to shield from cosmogenic muons — the primary non-antineutrino background source. Geoneutrino signals are often reported in terrestrial neutrino units (TNU), with one TNU equal to one IBD detection per 1032 free protons ($\approx$ one kiloton of liquid scintillator) in one year with 100% detection efficiency of a geoneutrino detector. This unit accounts for differences in detector size and efficiency.

Antineutrino detectors are presently counting geoneutrinos at KamLAND (Kamioka, Japan; 1 kton) and Borexino (Gran Sasso, Italy; 0.3 kton); future detectors include SNO+ (1 kton; online 2019) in Sudbury, Canada, JUNO (20 kton; online 2022) in Guangdong, China, and Jinping (4 kton; unknown start date) in Sichuan, China. Although the Earth is emitting $\sim 10^{6}$ $\bar{\nu_{e}}/$cm${}^{2}/$s, few events are detected annually at KamLAND (∼14/yr) and Borexino (∼4/yr) because of the combined effects of the inverse square law (intensity $\propto$ 1/distance2) and the neutrino’s small interaction cross section. Of these detected events, \citeAaraki2005 estimated that 50% are derived from U and Th in the HPE-enriched upper continental crust within $\sim 500$ km of the detector (known as “near-field crust”). The remaining signal comes equally from the rest of the continental crust (i.e. the “far-field crust”; $>500$ km) and the mantle. Consequently, the signal from the crust ($\sim 40$ km thick) overpowers that from the more massive mantle ($\sim 2900$ km thick). In order to determine the mantle contribution to the geoneutrino flux — and therefore the amount of U and Th within the mantle — it is necessary to have precise and accurate estimates of the crustal contribution. Current measurements from KamLAND and Borexino do not yet allow for discernment between different BSE models (Figure 1).

Endeavors to understand the crustal signal constitute a major effort in geoneutrino research. Models of the global crust employ a combination of global and regional seismic and gravity data with extrapolation to areas with minimal available data. Similarly, geochemical data from global compilations representative of sediment, upper, middle, and lower crustal layers are combined with these physical models. These joint models make simplifying assumptions regarding the bulk composition of each layer of the continental crust. Only recently did these geochemical model estimates include lateral and vertical spatial variability, although not in all layers of the crust <e.g.,¿huang2013. A further confounding issue is that there exist different geophysical models of the bulk lithosphere and each yields slightly different geoneutrino estimates. From the variance between these models and lack of direct measurement, it is unclear which model provides a more accurate representation of the crust. For example, estimated signals at KamLAND can vary by as much as 40%, with particular distinction between models produced respectively by the geoscience and physics communities (Figure 2).

This paper aims to update the geoneutrino model of \citeAhuang2013 by application of their geochemical methods to three different geophysical reference models, including a more recent physical model. First, the geophysical structure and properties of the crust from three different geophysical models are described. Geochemical methods are applied to these geophysical models to yield a 3D description of the amount and distribution of [U,Th,K] in the crust. Finally, these combined geophysical and geochemical models are used to calculate the geoneutrino signal at six detector locations. Uncertainties for all input parameters are defined, correlated, and propagated using Monte Carlo methods. Finally, geochemical and geoneutrino signal modeling results are compared and discussed.

2 Geophysical Model

The lithosphere is the outer rigid silicate shell of the Earth and is composed of the thick ultramafic lithospheric mantle and its overlying, relatively thinner and more felsic crust. This crust is either mafic oceanic or intermediate continental crust. In the past few decades authors moved from characterizing the vertical profile of the continental crust as a single bulk layer to a combination of an upper crust and lower crust <e.g.,¿hacker2015, an upper, middle, and lower crust <e.g., CRUST5.1;¿mooney1998, or as a continuously changing density structure <e.g., GEMMA1.0;¿reguzzoni2015. Furthermore, authors have begun to laterally characterize the crust as an assemblage of different groups with similar tectonic and seismic structure, something that has previously been performed <e.g.,¿pakiser1966 but rarely with regard to both lateral and vertical variations. Modern global physical models are constructed using global and local seismic studies, gravity surveys, or a combination of both. The primary utility of global geophysical models is for crustal corrections in mantle tomography Mooney \BOthers. (\APACyear1998) and interpretation of gravity measurements <e.g.,¿reguzzoni2015. These models also serve as the basis for geochemical and geoneutrino modeling.

CRUST5.1 Mooney \BOthers. (\APACyear1998) provided seismic compressional ($V_{P}$) and shear wave ($V_{S}$) velocity averages for regions with similar crustal structure in order to extrapolate to regions where there was limited or no data Mooney \BOthers. (\APACyear1998). This model defined the crust as three layers (upper, middle, and lower) with a sediment layer on top. The subsequent family of geophysical models built on CRUST5.1’s data and methods have a generational increase in input data, resulting in resolution changes from CRUST5.1 (5°lat $\times$ 5°lon, i.e., 2592 tiles laterally) $\rightarrow$ CRUST2.0 <2° $\times$ 2°, 16200 tiles;¿bassin2000 $\rightarrow$ CRUST1.0 <1° $\times$ 1°, 64800 tiles;¿laske2013). Of note, CRUST2.0 and CRUST1.0 provide limited transparency of the model inputs or methodology aside from AGU and EGU abstracts, respectively. The most recent iteration of the ”CRUST” family of models, LITHO1.0 ($\sim$1° $\times$ 1°) Pasyanos \BOthers. (\APACyear2014), perturbed CRUST1.0 parameters (density, seismic speeds $V_{P}$ and $V_{S}$, layer thickness) to fit a global surface wave dataset and for the first time included lithosphere-asthenosphere boundary (LAB) depths. LITHO1.0 reports $V_{P}$ and $V_{S}$ as a continuous probability distribution, while CRUST2.0 and CRUST1.0 report a few discrete velocities within each layer. Recently, \citeAolugboji2017 showed that LITHO1.0 displays a smaller misfit to high-resolution seismic studies in the United States compared to CRUST1.0. CRUST2.0, CRUST1.0, and LITHO1.0 are adopted as the physical basis for our geoneutrino modeling in order to test their effect on modeled outputs. Model parameters relevant to this study are the layer thickness, density, and $V_{P}$. These geophysical models are not combined with the more recent GEMMA1.0 (as was done by \citeAhuang2013) due to ambiguity regarding GEMMA1.0’s dependence on the adopted density structure of the crust. However, this also means that we do not have a measure of uncertainty on crustal thickness and must rely on estimates from previous studies (see Section 4). A comparison of the mass and thickness of layers within these models is shown in Table 1, Figure 3 (density and mass comparisons), and Figure 4 (thickness comparisons).

Previous geoneutrino reference models were generally built upon the most recent model available. The first model post the start of geoneutrino measurements at KamLAND was that of \citeAmantovani2004 built using the Preliminary Reference Earth Model <PREM;¿dziewonski1981, followed by the models of \citeAenomoto2007 and \citeAfiorentini2012 for the geoneutrino signal at KamLAND and Borexino (built on the geophysical model CRUST2.0). \citeAhuang2013 is built on a combination of the physical models CRUST2.0, CUB2.0 Shapiro \BBA Ritzwoller (\APACyear2002), and GEMMA1.0 and has been used in the prediction of signals at SNO+ Huang \BOthers. (\APACyear2014); Strati \BOthers. (\APACyear2017) and JUNO Strati \BOthers. (\APACyear2015), among others. Finally, \citeAsramek2016 calculated the signal at Jinping with a model built upon CRUST1.0 Laske \BOthers. (\APACyear2013).

For the purposes of this study LITHO1.0 is re-meshed from its original icosahedron equal-area tesselation (level 7; with $\sim$1° spacing) onto a 1°lat $\times$ 1°lon grid. Similarly, CRUST2.0 is re-meshed into a 1° $\times$ 1° grid from its previous 2° $\times$ 2° grid. Re-meshing was necessary to more easily compare modeled outputs. We split the LITHO1.0 and CRUST1.0 tiles into continental and oceanic crust using the crust type characterization from CRUST1.0 and split CRUST2.0 by its own classification. For CRUST1.0 and CRUST2.0, oceanic crust is defined as types A0, A1, B-, V1, Y3 or all A and B types, respectively; these correspond to crustal types labeled ”oceanic”. The final models contain sediment layers (two for CRUST2.0 and 3 for CRUST1.0 and LITHO1.0), upper, middle, and lower crust, and lithospheric mantle. Sedimentary thickness within the ”CRUST” family are from a digitization of energy industry datasets Laske \BBA Masters (\APACyear1997).

Continental lithospheric mantle LAB depths are provided by LITHO1.0; for CRUST2.0 and CRUST1.0 the LAB is set to 175 $\pm$ 75 km depth following \citeAhuang2013. Although this study focuses on the lithosphere and continental crust, we include a convecting mantle in order to provide an estimate of the total geoneutrino signal expected at a detector. For the convecting mantle we define the lower 750 km as the Enriched Mantle, which equals $\sim$19% of the mantle by mass Dziewonski \BBA Anderson (\APACyear1981); Arevalo \BOthers. (\APACyear2013). PREM densities are adopted for the mantle Dziewonski \BBA Anderson (\APACyear1981). The mass of the BSE reported in Table 1 is a combination of PREM for the convecting mantle and either CRUST2.0, CRUST1.0, or LITHO1.0 for the lithosphere.

[FIGURE:]

3 Geochemical Model

The Earth’s crust is composed of mafic oceanic crust (on average $\sim 7$ km thick) and relatively more felsic continental crust (avg. $\sim 35$ km thick). The geochemistry of the continental crust has been explored extensively over the past century, particularly the accessible upper crust. Generally the crust is thought to decrease in $SiO_{2}$ with depth from 67% in the upper crust to 53% in the lower crust Rudnick \BBA Gao (\APACyear2014). These observations are from surface exposures of the upper, middle, and lower crust as well as xenolith samples. Furthermore, it has long been observed that increasing metamorphic grade cannot account for the increase in the observed $V_{P}$ with depth, rather requiring a compositional increase in maficity Christensen \BBA Mooney (\APACyear1995). The relationship between seismic velocity and chemical composition has been observed empirically, due to the chemical effect on the shear and bulk modulus of which $V_{P}$ relies Christensen (\APACyear1965); Pakiser \BBA Robinson (\APACyear1966); Christensen \BBA Mooney (\APACyear1995). However, observed $V_{P}$ of different lithologies are non-unique, requiring further assumption to have any meaningful result when inverting $V_{P}$ for composition. The methodology of \citeAhuang2013 is adopted here, which assumes the continental middle crust is composed of amphibolite and the continental lower crust of granulite metamorphic rocks. These assumptions are based on surface exposed crustal cross-sections and xenolith data Rudnick \BBA Gao (\APACyear2014). \citeAhuang2013 correlated $V_{P}$ and $SiO_{2}$ and $SiO_{2}$ and [U, Th, K] abundance of amphibolite and granulite by a linear relationship between mafic and felsic endmembers. This correlation is characterized by a simple mass balance defined by the following equations:

[TABLE]

[TABLE]

[TABLE]

where $V_{f}$, $V_{m}$, $a_{f}$, and $a_{m}$ are the $V_{P}$ and abundance of the felsic and mafic end-members of amphibolite (middle crust) and granulite (lower crust), $V_{model}$ is the $V_{P}$ provided by the geophysical model, and f and m are the mass proportions of felsic and mafic endmembers. The velocities of the amphibolite and granulite datasets are derived from laboratory measurements at room temperature and 600 MPa Huang \BOthers. (\APACyear2013). These measurements are temperature corrected ($-4\times 10^{-4}$ km/s/°C) following a mid-range geotherm ($q_{s}=60$ mW/m2) Pollack \BBA Chapman (\APACyear1977); Turcotte \BBA Schubert (\APACyear2014)) and pressure corrected ($2\times 10^{-4}$ km/s/MPa) from pressures calculated from the parameters provided by the geophysical model ($\mathrm{d}P=\rho g\mathrm{d}h$) Huang \BOthers. (\APACyear2013). Temperature and pressure corrections are from empirical studies Christensen \BBA Mooney (\APACyear1995); Rudnick \BBA Fountain (\APACyear1995). In cases where the $V_{P}$ from the geophysical model is larger than the mafic endmember or smaller than the felsic endmember, we set $f$ to be 0 or 1, respectively. In LITHO1.0, this occurs in $\sim$11% (36% $f\!=\!0$ and 84% $f\!=\!1$, proportionally) of middle crustal and $\sim$42% (98% $f\!=\!0$ and 2% $f\!=\!1$, proportionally) of lower crustal continental tiles. Model $V_{P}$ outside the range of endmembers also occurs in CRUST2.0 and CRUST1.0, with the middle (lower) crust having 11% (23%) and 0% (64%), respectively. This method is especially problematic for the lower crust, but for the purposes of this study the impact has a negligible effect on the geoneutrino signal as the lower crust is generally depleted in [U,Th,K] and is farther from the geoneutrino detectors than the upper or middle crust (see Section 5).

HPE abundance for the upper crust is not calculated from a $V_{P}$ and composition correlation because near surface processes strongly affect observed $V_{P}$ and because we cannot simplify the upper crust as a single metamorphic grade (as we do for the middle and lower crust). Without the simplifying assumption of a single metamorphic grade there are too many non-unique lithologies that could be present. Instead values calculated from a meta-analysis of previously published upper crust estimates is adopted Rudnick \BBA Gao (\APACyear2014), which equates to a sigma-mean value rather than a standard deviation. Although this uncertainty is not consistent with our overall adoption of 1-sigma uncertainties in this study, the analysis of \citeArudnick2014 currently provides the best estimate of upper crust composition.

Sediment compositions are from the GLOSSII model of subducted sediments, as these material would be the remnants of continental weathering and marine activity Plank (\APACyear2014). This assumption is acceptable as the location of sedimentary basins in the geophysical model is generally along coastlines. The continental lithospheric mantle is characterized by a suite of xenoliths compiled by \citeAhuang2013. For the oceanic crust we adopt the bulk oceanic crust composition from \citeAwhite2014. The lithospheric mantle below oceanic crust is assumed to have the same composition as the depleted mantle.

The BSE composition is calculated from the U abundance given by \citeAmcdonough1995, whereas Th and K are calculated from observed Th/U and K/U ratios Arevalo \BOthers. (\APACyear2009); Wipperfurth \BOthers. (\APACyear2018). Similarly, a conservative upper mantle (LAB to 750 km above the core-mantle boundary) U abundance is assumed based on the compilation from \citeAarevalo2013 with the Th and K abundance calculated from observed Th/U and K/U values from mid-ocean ridge basalt samples Arevalo \BOthers. (\APACyear2009); Wipperfurth \BOthers. (\APACyear2018). The abundances within the enriched mantle are the remainder from subtraction of the lithosphere and upper mantle from the BSE composition (Table 2).

4 Uncertainties and Correlation

The attribution of uncertainty is the most difficult and time consuming part of numerical modeling. When available, uncertainties directly reported by the data source are adopted. When not available, relative uncertainties in agreement with literature estimates are adopted and correspond to 1$\sigma$, except for abundances in the upper crust. [U,Th,K] abundances display log-normal distributions, a trait shared with other incompatible elements Ahrens (\APACyear1954); McDonough (\APACyear1990); Wipperfurth \BOthers. (\APACyear2018). Model outputs, which depend on asymmetrical abundance distributions, are reported as geometric mean $\pm$ 1$\sigma$ uncertainties. 5% uncertainty is ascribed to $V_{P}$, which is primarily based on the recent comparison of a high-resolution surface wave model of the US with CRUST1.0 and LITHO1.0 Olugboji \BOthers. (\APACyear2017), although this is in agreement with earlier estimates Mooney \BOthers. (\APACyear1998). Because density is derived from $V_{P}$ in the ”CRUST” family of models, 5% uncertainty is also applied to density. 12% uncertainty is ascribed to crustal thickness based on the comparison of CRUST2.0, CUB2.0, and GEMMA1.0 by \citeAhuang2013. Uncertainty on mantle thickness are not incorporated in order to decrease computation time, although the effect on geoneutrino signal is expected to be negligible as the uncertainty on the mantle mass is $\sim 0.5\%$.

The correlation of uncertainties is equally as important as absolute uncertainty attribution and largely dictates the relative magnitude of the uncertainty on the modeling output. For example, no correlation in thickness of layers, and thereby no correlation of the mass, results in negligible uncertainty on the estimated mass of the crust ($\ll$1%). This result is unrealistic as we do not know the crustal mass to such precision. Thickness is correlated vertically across all layers and laterally across all 1° $\times$ 1° tiles; in this way the crust is modeled when it is most massive (all layers as thick as possible) and least massive (all layers as thin as possible). A partial uncertainty would be ideal but is impossible without in-depth knowledge of the correlation of input data used to construct the geophysical model (e.g., the foundational data of LITHO1.0 and their correlations). Neither $V_{P}$ nor density are correlated between the middle and lower crust as it is unclear how $V_{P}$ and density would shift given a change in one of the layers. Similar to thickness, $V_{P}$ and density are laterally correlated within a layer. Abundances are correlated within each voxel as there is a longstanding observation of a general correlation between [U,Th,K] due to their similar incompatibilities during mantle and crustal melting. This correlation also applies to the endmember abundances for amphibolite and granulite. However, abundances are not correlated vertically between layers because the samples used to estimate the abundances of the upper, middle, and lower crust are not the same samples and therefore variability within each dataset is independent of the others. There is, however, an inherent correlation between the middle and lower crust in our model because of the correlation of $V_{P}$ within these layers from the geophysical model coupled with our geochemical methods. Abundances are correlated laterally within a layer, which coupled with lateral correlations of $V_{P}$, density, and thickness results in models when layers have the most massive and least massive amount of HPE possible. All correlations described above assume 100% correlation between parameters. Our correlation approach is conservative; removal of lateral correlations would significantly decrease geoneutrino signal uncertainties but is not obviously justified at this time.

5 Numerical Model and Geoneutrino Signal

The geophysical and geochemical information are combined into a coherent model, with attributes assigned for each 1° $\times$ 1° cell for each layer. A Monte Carlo simulation is used to propagate uncertainty with 3$\times$104 iterations, as at this level we reach stability of outputs. In each iteration the input variable’s probability density function (PDF) is randomly sampled. If the parameter is correlated with another, then both parameters are sampled from the same part of the PDF. This method allows for easy combination of normal and log-normal distributions. The code is written for parallel computing in MATLAB using either LITHO1.0, CRUST1.0, or CRUST2.0, taking $\sim$8 hours on a 4 core (8 thread) Intel CPU for $3\times 10^{4}$ iterations of the Monte Carlo. In this way we have provided an environment to use and test the geoneutrino response given different geophysical models within a self-consistent framework. The effect of the number of iterations on the final signal was tested by performing 10 repeated calculations with $20\times 10^{4}$ iterations each, with relative variation in the central value and 1-sigma uncertainties of 0.1% and 0.5%, respectively.

For each voxel (3D pixel) within a layer, the mass was calculated from the the density and thicknesses provide by the geophysical model. Multiplication of the mass by HPE abundance yields the mass of HPE within each voxel. The heat production is calculated using updated decay and heat output parameters from \citeAruedas2017. Finally, the geoneutrino signal can be calculated for U and Th. Iterations where the enriched mantle is more depleted in HPE abundances than the depleted mantle, which occurs in $\sim$25% of iterations, are discarded. This avoids conflicts with observed OIB samples with enriched HPE abundances compared to MORB samples Arevalo \BOthers. (\APACyear2013) as well as avoiding an enriched mantle with no HPE.

5.1 Calculating Signal

We calculate the predicted geoneutrino signal at a specified detector location defined by latitude, longitude, and depth below surface (see Supplementary Information). The energy spectrum of detected antineutrinos is written <e.g.,¿dye2012 as

[TABLE]

where the meaning of symbols and their units are explained in Table 3.

This equation accounts for the number of decays of each isotope, the emitted neutrino energy spectrum, the distance the isotope is from the detector, and the probability of a neutrino oscillation over that distance. The oscillation parameter $P_{ee}$ is calculated from the most recent values from Capozzi \BOthers. (\APACyear2017). The signal is calculated and propagated for every 75 keV across the IBD-detectable geoneutrino emission spectrum ranging from $\sim 1800$ keV to $\sim 3300$ keV, as calculated by \citeAenomoto:nuspec. The output from equation 4 is in TNU.

As the geoneutrino signal is dependent on the distance from the detector in inverse square ($|\vec{r}-\vec{r}\prime|^{2}$) it is necessary during modeling to increase the mesh resolution for tiles near the detector, otherwise the signal would be calculated from only the center of each 1° $\times$ 1° voxel. The mesh resolution is increased until the final signal reaches stability. This results in modeled cell sizes increasing from $\sim$ 100 m3 next to the detector, 10 km3 at 500 km distance, and the initial cell size after 1000 km. The re-meshing is a necessary step in the calculation of the geoneutrino signal but is computationally intensive, increasing computation time by $\sim$40%.

6 Results and Discussion

This study provides an updated reference model to \citeAhuang2013 for the abundance of heat producing elements in the lithosphere, their heat production, the geoneutrino signal at current and future detectors, and a comparison of these outputs from three commonly used geophysical models (Table 4). The crustal abundances and geoneutrino signal are a result of the complex amalgamation of laterally and vertically variable seismic velocity, layer thickness, layer density, and [U,Th,K] abundance. These parameters change from cell to cell in a convoluted way, making analysis of the specific cause of differences in the geoneutrino signal and crustal abundances between geophysical models difficult. For these reasons the following section will attempt to convey some of the major differences when using different geophysical models. A broad overview of different lithospheric thickness models is provided by \citeAsteinberger2016, which complements this study.

6.1 Geoneutrino Signal

The geoneutrino signal results at six detectors are reported in Table 4, which includes signal from the bulk continental crust and total expected signal from the BSE. A full breakdown of each signal is included in the Supporting Information. This study finds that the near-field contributes $\sim$40% of the signal at the detector, with the remaining signal from the rest of the crust ($\sim$35%) and the mantle ($\sim$25%). To illustrate the similarity between signals when using CRUST2.0, CRUST1.0, and LITHO1.0, we include in Table 4 the maximum difference in the central value as well as the Overlapping Coefficient (OVL), which measures the degree with which the distributions overlap Inman \BBA Bradley Jr. (\APACyear1989). In our case the OVL is calculated from

[TABLE]

where $f_{M_{1}}$ and $f_{M_{1}}$ are normalized probability density functions for geoneutrino signal $S$ calculated with Earth models $M_{1}$ and $M_{2}$, and $M_{i}$ is selected from CRUST2.0, CRUST1.0, or LITHO1.0. The OVL can only reach values from 0 (no overlap) to 1 (complete overlap). The results of this study show that the choice of global physical model has only a small effect on the predicted crustal geoneutrino signal (Table 4). Estimates for the crustal geoneutrino signal at KamLAND, SNO+, and Jinping show larger variation between predictions than the other models (15% maximum variation in the central value). OVL values are always greater than 70% overlap, with average overlap greater than 90%. The OVL values suggests more similarity between LITHO1.0 and CRUST1.0 compared to CRUST2.0. This is in contrast to the apparent similarity in regional masses between CRUST2.0 and CRUST1.0 observed in Figure 3.

[FIGURE:]

The measured signal at Borexino <$43.5^{+11.8}_{-10.4}$ TNU;¿agostini2015 and KamLAND <$34.9^{+6.0}_{-5.4}$ TNU;¿[]watanabe2016 are plotted against the predicted lithospheric signal from LITHO1.0 in Figure 5 in order to interpret the observed signals for BSE radiogenic power. Figure 5 represents a two-component system whereby the total measured signal (from the BSE) is comprised of the lithospheric and mantle signals. The slope is therefore defined equal to unity (i.e., y = x + b). The y-intercept, which represents the signal from the mantle, is calculated using Monte Carlo methods in order to account for asymmetric X and Y uncertainties. The predicted lithospheric signal for KamLAND and Borexino are correlated during the y-intercept Monte Carlo as most parameters are correlated within the geoneutrino reference model. Conversion of mantle TNU signal to BSE radiogenic power uses the relationship derived using the geoneutrino prediction model (similar to Figure 1 which shows total geoneutrino signal vs. radiogenic power in BSE). The estimated BSE radiogenic heat production of $20.3\pm 12.2$ TW is consistent within 1-sigma of the low, medium, and high-H BSE models with the central value coincident with the medium-H model ($20\pm 4$ TW). Furthermore, the prediction given here using LITHO1.0 is consistent with that predicted by the \citeAsramek2016 model ($21\pm 10$ TW) when using the updated \citeAwatanabe2016 KamLAND measurement (personal communication). Based on Figure 5, an observed mantle signal of 8.5 TNU with uncertainty below some $\pm 6$ TNU is needed to discriminate between the different BSE compositional models. This uncertainty corresponds to a 40% reduction in the measured total and predicted lithospheric signal uncertainty at KamLAND and Borexino.

6.2 HPE Abundance and Heat Production

Abundances described in Table 2 for all layers are comparable (within 1$\sigma$) to both \citeArudnick2014 and \citeAhuang2013. However, uncertainties reported in the middle and lower crust are larger than those reported by \citeAhuang2013, likely due to differences in how these values were correlated. K/U and Th/U of the layers of the crust and bulk crust are also in agreement at 1$\sigma$ with previous estimates Huang \BOthers. (\APACyear2013); Rudnick \BBA Gao (\APACyear2014); Wipperfurth \BOthers. (\APACyear2018).

The bulk-continental crust heat production is negligibly different between the geophysical models (Table 2) and is consistent with previous estimates Mareschal \BBA Jaupart (\APACyear2013); Huang \BOthers. (\APACyear2014); Rudnick \BBA Gao (\APACyear2014). Furthermore, an estimate of bulk continental heat production of 7.0${}^{+2.1}_{-1.6}$ TW (using LITHO1.0) encompasses average heat production in stable Precambrian (5.9 $\pm$0.6 TW) and Phanerozoic (8.3 $\pm$1.0 TW) crust estimated from heat flow studies Mareschal \BBA Jaupart (\APACyear2013). The similarity of continental crust heat production when using CRUST2.0 and LITHO1.0 is in contrast with the crustal signal similarity of CRUST1.0 and LITHO1.0 (Table 4) or the apparent crustal mass similarity of CRUST2.0 and CRUST1.0 observed in Figure 3 panel C. The similarity in crustal heat production of CRUST2.0 and LITHO1.0 is simply a reflection of their similar bulk crustal mass compared to the lower mass of CRUST1.0 (Table 1).

6.3 Uncertainties

The reported uncertainties for crustal abundances of the HPE in Table 2 are larger than that reported in \citeAhuang2013, primarily as a result of increased correlation of abundances and geophysical uncertainties. Increased correlation causes an inverse effect on uncertainty for the bulk crustal K/U and Th/U, resulting in less uncertainty than described by \citeAhuang2013. Some previous geoneutrino signal predictions did not include uncertainty on geophysical inputs, including crustal thickness <e.g.,¿enomoto2007,fiorentini2012,sramek2016, meaning these studies report smaller uncertainty on the geoneutrino signal. Similar to signal, the uncertainty on continental crustal heat production is $\sim$50% larger than that estimated by \citeAhuang2013 ($6.8^{+1.4}_{-1.1}$ TW). In general, geophysical uncertainties (on thickness, density, and $V_{P}$) account for $\sim 30\%$ of the signal uncertainty, with the remaining proportion ($\sim 70\%$) from the geochemical inputs. When only considering the geophysical uncertainty, as the same geochemical method was applied to each geophysical model, the calculated crustal signals are still negligibly different.

7 Future Developments

The results from this study highlight the convergence of our understanding of the physical nature of the continental crust and the negligible effect of geophysical model on the predicted geoneutrino signal. Due to the strong effect of geochemical variability and distance from the detector on the predicted geoneutrino signal, future geoneutrino modeling needs to focus on better understanding of geochemical variability and a focus on near-field modeling. Below we highlight a few areas which require further exploration in order to better understand the crustal component of the geoneutrino signal.

7.1 Upper Crust Geochemical Uncertainty

The attribution of uncertainty on the concentration of heat producing elements in the upper crust (sigma-mean) in most geoneutrino studies (including this one) is not consistent with uncertainty on other parameters (1-sigma), including abundances in other layers. Often cited studies on the upper crust composition report sigma-mean Rudnick \BBA Gao (\APACyear2014) or Median Absolute Deviation <MAD;¿[]gaschnig2016, both of which report smaller error estimates than the standard deviation. Any regional or global geoneutrino modeling should adopt consistent error estimators, be it standard deviation or MAD, or justify the treatment of the upper crust separately from other portions of the crust. Upper crust U abundance used in this study is 2.7 $\pm$ 0.6 $\mu g/g$ (sigma-mean) from \citeArudnick2014, similar to that from \citeAgaschnig2016 of 2.66 $\pm$ 0.87 $\mu g/g$. The standard deviation of the glacial diamictite data (assumed to sample large portions of the upper crust) from \citeAgaschnig2016 is 150% of the MAD. Because the upper crust is the dominant heat and geoneutrino emitter, consistent uncertainty estimates on abundances could significantly increase uncertainties on the predicted crustal geoneutrino signal. For example, attribution of abundances in the upper crust using 1-sigma values from a log-normal fit of \citeAgaschnig2016 data yields a total signal at KamLAND of 36.4 ${}^{+9.8}_{-7.7}$ TNU (compared to 37.9 ${}^{+6.6}_{-5.6}$ TNU reported in Table 4) using LITHO1.0. Use of 1-sigma uncertainty for the upper crust abundances therefore results in an increase of uncertainty of $\sim$40% in the predicted signal, an observation consistent with other studies <e.g.,¿takeuchi2019.

7.2 $V_{P}$ and Density under Himilayas and Andes

Future studies with the goal to update LITHO1.0 should aim to better understand the density and seismic structure beneath the Himalayas and Andes regions. The ’CRUST’ family of models predicts felsic-like densities in the entire crustal column in these regions (including the expected mafic lower crust; Figure 3). This prediction is in conflict with some regional studies <e.g.,¿monsalve2008,bai2013 although in agreement with other global studies <e.g.,¿hacker2015 and regional seismic models Agius \BBA Lebedev (\APACyear2017); Gilligan \BBA Priestley (\APACyear2018). Additionally, \citeApasyanos2014 perturbed CRUST1.0 parameters up to 5% to fit observed surface waves to make LITHO1.0, which a comparison of CRUST1.0 and LITHO1.0 $V_{P}$ in the Himilayan region shows that this perturbation reached saturation (i.e. 5% change between LITHO1.0 and CRUST1.0). If the model was not limited to 5% change the output would be more felsic than current. Because of the saturation of the method used by \citeApasyanos2014, LITHO1.0 in its current state may not be in agreement with observed regional data and should be revisited. Similar phenomena are observed under the Andes Mountains in South America <e.g.,¿lucassen2001. The Himilayas are particularly relevant for the geoneutrino prediction signal at Jinping and to a lesser extent, JUNO, while the Andes will be relevant to the future ANDES underground laboratory Bertou (\APACyear2012).

7.3 Geochemical Method for Middle and Lower Crust

The geochemical method of conversion of $V_{P}$ to [U,Th,K] does not work well in the lower crust (and to some degree in the middle crust) as the model often surpasses the endmember condition. This is a problem due to simplifying assumptions (i.e. assuming endmembers) indicating a need for a more sophisticated modeling space. A bivariate probability analysis of the available amphibolite or granulite samples with measured $V_{P}$ could avoid this problem (see Supporting Information for details). This analysis would not assume any specific relationship between $V_{P}$ and $SiO_{2}$ (unlike the linear relationship we assumed in this study) but instead would use the bivariate probability of $V_{P}$ and $SiO_{2}$ from the dataset. Currently there are too few samples ($\sim$100–150) with measured $V_{P}$ to create a robust analysis. Incorporation of a thermodynamic modeling software <such as Perple_X;¿connolly2005) would allow for the calculation of seismic wavespeed for $>$500 samples, which would significantly increase the robustness of a bivariate probability analysis. Calculated $V_{P}$ using Perple_X for samples used in \citeAhuang2013 by \citeAhacker2015 closely resemble laboratory measurements, indicating the viability of thermodynamically calculated $V_{P}$.

7.4 Near-field (regional) Modeling

Underestimation of the upper crust uncertainty or problems associated with middle/lower crust abundance calculations have less of an effect on predicted signals if a high-resolution regional geoneutrino model is combined with the global model <e.g.,¿enomoto2007,coltorti2011,huang2014,strati2017. Because of the distance dependence of the geoneutrino signal (see eq. 4) the regional area provides $\sim$40% of the geoneutrino signal at a detector as calculated in this study. High-resolution seismic and geochemical studies of the near-field therefore provide the most robust estimate of the signal at any detector location. This is exemplified in the studies of \citeAhuang2014 and \citeAstrati2017, who calculated larger uncertainties on the geoneutrino signal at SNO+ when they included a regional model compared to only using a global model. Although the uncertainty is larger than that estimated from the global model, their estimate includes more local information and is therefore likely to be a more accurate estimate of uncertainty.

Finally, evolution of our understanding of the physical and chemical nature of the bulk-crust is largely plateauing. This study showed that geoneutrino estimates using LITHO1.0 (2015), CRUST1.0 (2013), and CRUST2.0 (2001) yielded largely the same geoneutrino signal. Furthermore, the geochemical description of the crust from this study is consistent with the estimates by \citeArudnick2014, among others, which is itself consistent with \citeArudnick2003. Because of the natural variability of [U, Th, K] in the crust the uncertainty on bulk abundances is unlikely to change significantly. Joint inversions of multiple datasets, including teleseismic $V_{P}$ and $V_{S}$, surface waves, gravity, surface heat flow, and topography/geoid have the most potential for significant advances in our understanding of the structure and composition of the crust (e.g., \citeAafonso2016,afonso2019). Regardless, improving our understanding of the physical and chemical nature of the bulk-crust will not have the same impact on reduction in the geoneutrino signal uncertainty as improvements in the near-field. As highlighted by \citeAstrati2017, better understanding of the physical structure with depth is the most difficult hurdle to overcome in near-field modeling.

8 Author Contributions

S.A.W. modeled the geoneutrino signal and wrote the text, with significant input on model creation, data analysis, and the text from O.Š. and W.F.M.

Acknowledgements.

Model results and original MATLAB code is provided in the Supporting Material. Support for this study was provided by the University of Maryland Graduate School ALL-S.T.A.R. Fellowship and NSF EAPSI program Award #1713230 (to S.A.W), NSF grant EAR1650365 (to W.F.M.), and the Czech Science Foundation grant GAČR 17-01464S (to O.Š.). There are no financial or interest conflicts with this work. We thank Fabio Mantovani, Bedřich Roskovec, and Steve Dye for insightful comments and discussion.