Strong Sequestration of Hydrogen Into the Earth's Core During Planetary Differentiation

We explore the partitioning behavior of hydrogen between coexisting metal and silicate melts at conditions of the magma ocean and the current core–mantle boundary with the help of density functional theory molecular dynamics. We perform simulations with the two‐phase and thermodynamic integration methods. We find that hydrogen is weakly siderophile at low pressure (20 GPa and 2,500 K), and becomes much more strongly so with pressure, suggesting that hydrogen is transported to the core in a significant amount during core segregation and is stable there. Based on our results, the core likely contains ~1 wt% H, assuming single‐stage formation and equilibration at 40 GPa. Our two‐phase simulations further suggest that silicon is entrained in the core‐forming metal, while magnesium remains in the silicate phase. This preferred incorporation of silicon hints at an explanation for the elevated Mg/Si ratio of the bulk silicate Earth relative to chondritic compositions.


Introduction
The Earth's early accretion events, in particular magma ocean formation, had a great impact on the chemical and thermal evolution of the Earth. Segregation of metal from silicate to form a core during the magma ocean stage has removed elements besides iron and nickel from the mantle (e.g., Righter, 2003). Light elements such as H, C, O, and Si enter the core as a consequence of this process, which is supported by the observed density (ρ) deficit and sound velocity profiles in the core (e.g., Birch, 1952;Poirier, 1994;Wagle & Steinle-Neumann, 2019).
As the most abundant element in the universe, hydrogen is one of the candidates. Although its solubility in iron at atmospheric pressure (P) is very limited, the affinity of hydrogen to iron increases significantly under compression (Antonov et al., 1998;Badding et al., 1991;Fukai, 1984) as the Gibbs energy (G) of gaseous hydrogen increases rapidly with P. Solution of hydrogen in iron has been confirmed by neutron diffraction measurements (Ikuta et al., 2019;Machida et al., 2014) at high P, with H dissolving interstitially in Fe, expanding the lattice in a systematic way. The affinity to iron (siderophile behavior) at high P suggests that hydrogen is a strong candidate for the light element in the core. This is reflected by decades of work on the Fe-H system at high P, including melting point depression and phase diagram (Okuchi, 1998;Sakamaki et al., 2009), hydrogen partitioning experiments between silicate and iron (Clesi et al., 2018;Malavergne et al., 2019;Okuchi, 1997), and by extension studies of water-iron reactions (Fukai, 1984;Ohtani et al., 2005;Yagi & Hishinuma, 1995;Yuan et al., 2018). Recent diamond anvil cell (DAC) experiments reveal that a series of iron hydrides, e.g., FeH (Hirao et al., 2004;Isaev et al., 2007;Meier et al., 2019;Narygina et al., 2011), and compounds with higher H contents (Pépin et al., 2014(Pépin et al., , 2017 are thermodynamically stable at conditions relevant to the deep Earth. Okuchi (1997) pioneered hydrogen partition experiments between molten iron and silicate at P to 7.5 GPa in the multi-anvil (MA) press to explore its distribution during core-mantle differentiation. Okuchi (1997) suggested that > 95% of the H 2 O present in a hydrous magma ocean may have reacted with Fe to form FeH x , and~60% of the ρ-deficit of the core can be accounted for by hydrogen as a light element. However, recent metal-silicate partitioning experiments by Clesi et al. (2018) and Malavergne et al. (2019) found hydrogen to be lithophile at high P and T, which shifted the debate of hydrogen in the Earth's core. Determining the hydrogen content of the metal in such experiments remains challenging as FeH x decomposes during sample recovery from high P and T. Also, graphite was used as sample capsules, and recent DAC experiments  suggest that C and H are mutually exclusive in molten iron. The presence of carbon may reduce the hydrogen storage capacity of iron, resulting in a very small equilibrium constant K D between molten iron and silicate melt (1) with x met H (x sil H ) the mole concentrations of H in the metallic (silicate) melt. Hydrogen metal-silicate partitioning data are only available at moderate P up to 21 GPa due to experimental challenges (Clesi et al., 2018;Malavergne et al., 2019;Okuchi, 1997). Here, we present results from molecular dynamics (MD) simulations using density functional theory (DFT) based methods on the interaction between hydrogen-bearing silicate and metallic liquids, determining hydrogen partitioning at 20-130 GPa and 2,500-4,000 K, which covers the whole range of P-T conditions of metal-silicate equilibration/segregation during core formation (Li & Agee, 1996;Siebert et al., 2012).

Density Functional Theory Molecular Dynamics Simulations
We conduct first-principles electronic structure calculations based on Kohn-Sham (KS) DFT to obtain total energy, Hellmann-Feynman stresses and forces with a plane wave approach and periodic boundary conditions using the Vienna Ab Initio Simulation Package (VASP) (Kresse & Furthmüller, 1996;Kresse & Hafner, 1993) with a planewave energy cutoff of 450 eV. We use the projector augmented wave implementation and atomic files (Kresse & Joubert, 1999) and the generalized gradient approximation to the exchange and correlation potential (Perdew et al., 1996). Electronic KS-DFT states are calculated at the Brillouin zone center only. Molecular dynamics (MD) simulations are performed in the canonical N-V-T ensemble, where the number of atoms (N) and the volume (V) of the cell are fixed, and T is controlled by a Nosé-Hoover thermostat (Hoover, 1985;Nosé, 1984). For the calculation of Helmholtz energy (F) using thermodynamic integration simulations (Section 2.3), the Nosé-Hoover thermostat is not stable when the system approaches the ideal gas limit with small distances between atoms, and we use the Langevin thermostat instead (Allen & Tildesley, 2017), following the approach of Dorner et al. (2018). After convergence of the electronic cycle, atoms are advanced with time steps of Δt = 0.5 fs due to the fast dynamics of hydrogen andas for the thermostatwe change this to smaller Δt = 0.2 fs as we approach ideal gas behavior for TI. Initial structures of the liquids are generated by performing MD calculations at T = 8,000 K, which is then changed to the target T. We check that the cells considered are indeed in the molten state by monitoring the mean-square displacement (MSD) of all species and partial radial distribution functions (RDF) between them.

Two-Phase Simulation
Two-phase simulations model element partitioning directly, mimicking an experiment by putting metallic and silicate phases in direct contact in a simulation cell (SI, Figure S1), and then tracking the distribution of hydrogen atoms between the two liquids as the system approaches equilibrium. The two-phase technique has been widely used to study the melting T of different materials at high P (e.g., Alfè, 2005Alfè, , 2009 Figures S1a, S1b); a small gap between the metallic and silicate surfaces has to be introduced to avoid that distances between some atoms near the interface become too small. Even though hydrogen is likely present in its reduced oxidation state in silicate liquids under reducing conditions of deep magmas (Hirschmann et al., 2012), we also explore partitioning of hydrogen under oxidizing conditions by using the composition Mg 37 Si 37 O 114 H 6 -Fe 150 O 3 H 6 (MSHO) with 353 atoms (SI, Figure S1c). Configurations from the DFT-MD trajectory for each of the conditions considered (SI ,  Table S1) are analyzed in terms of the local atomic structure and spatial distribution of all atoms, in particular, their number density profile in the simulation box along the axis perpendicular to the phase boundary.
Although two-phase simulation are very straightforward, there are significant limitations to this method which highlight the need for additional calculations: i. It is not appropriate for studying the partitioning behavior of atoms with similar or larger masses than those of the main elements constituting the silicate and metallic liquids (i.e., Mg, Si, O, and Fe).
Hydrogen diffuses at least one order of magnitude faster than other kinds of atom in the melts (e.g., Posner & Steinle-Neumann, 2019). ii. To observe partitioning behavior, simulation T have to be moderate (here, 2,500-3,500 K), which ensures that hydrogen equilibrates between two coexisting phases, but prohibits extensive mixing of silicate and metallic phases. The mobility of heavy atoms can be enhanced substantially at high T, which leads to complete mixing of silicate (or silica) and metallic liquids, as observed in previous DFT-MD simulations in the Fe-Si-O ternary (Huang et al., 2019). iii. The two-phase simulation method can only provide qualitative results on metal-silicate partitioning of hydrogen (i.e., determine whether it is lithophile or siderophile) without the possibility to compute partitioning coefficients.

Thermodynamic Integration
We use the following reaction to evaluate hydrogen partitioning between the metallic liquid and silicate liquid: where we exchange H 2 O rather than H 2 for computational convenience and as justified by the similar H behavior of the two-phase simulations for the MSH and MSHO systems (Section 3.1), and no inference is made concerning the nature of H speciation. The Gibbs energy change (Δ r G) of this reaction is The equilibrium constant of hydrogen defined in terms of the mole concentrations of H in the metallic and silicate melts (eq. 1) can then be calculated as where k B is the Boltzmann constant.
DFT-MD simulations give potential energy (U) and pressure (P) but not entropy (S) of the system; combining DFT-MD with the thermodynamic integration method provides a general scheme for computing the Helmholtz energy (F) difference ΔF = F B -F A of two systems containing the same number of atoms which are governed by different potential energy functions U B (target system) and U A (reference system) (e.g., Alfè et al., 2000;Vočadlo et al., 2008;Wahl & Militzer, 2015;Xiong et al., 2018). For the metal and silicate liquids considered in this study, we use the ideal gas as reference, for which F can be calculated analytically. We perform a direct thermodynamic integration from ideal gas to the fully DFT interacting system at fixed V and T (Dorner et al., 2018;Rang & Kresse, 2019;: where the notation ⟨A⟩ λ signifies the thermodynamic average of quantity A over an ensemble generated by the Hamiltonian H(λ) for which the potential energy U(λ) term is a linear mixture of U A and U B (U(λ) = (1 -λ) U A + λ U B ), with λ a coupling constant that connects the reference to the target system. The advantage of this method is that it expresses ΔF in terms of an ensemble average, whichunlike F or Gcan be calculated directly in a simulation (Frenkel & Smit, 2001): i. A system with potential energy function U(λ) is simulated at a series of λ values. The derivative of U with respect to λ is then averaged over thermodynamically equilibrated configurations. ii. Separate standard DFT-MD simulations are performed to locally determine isothermal EoS of the system (Section 2.1; SI, Figure S2). iii. Numerical integration with respect to λ is carried out to obtain ΔF between the ideal gas and KS-DFT system. iv. Once ΔF and absolute F are computed by adding the ideal gas contribution F ig , G can be calculated by adding a PV term, i.e., G = F + PV. The EoS obtained in step (ii) are used to calibrate all the calculations to constant P, following the fundamental thermodynamic relation for G along an isotherm: DFT-MD simulations with various λ parameters for each of the compositions indicated in reaction 2 are performed at 20 GPa (2,500 K), and 40 GPa and 130 GPa (4,000 K) (SI, Figure S3) as indicated in Tables S2-S5 in the SI.

Two-Phase Simulations
Representative snapshots of the two-phase simulation results for the MSH, MSHC and MSHO compositions at two P-T conditions are shown in Figure 1 (SI , Table S1), animations for the development of the cells at all P-T conditions can be found in the supporting online material (SI, Movie S1). Simulations start from initial configurations where six hydrogen atoms each are located in the silicate (H S ) and metallic (H M ) liquids. At all P-T conditions considered, we observe the diffusion of hydrogen atoms from the silicate to the metallic phase, leading to an accumulation of hydrogen in iron (SI, Figure S4 and Movies S2-S4) that is also documented by the time-development of the partial radial distribution function (RDF) g(r) between both hydrogen species (H S and H M ) and Fe and O, i.e., g OHS r ð Þ, g OHM r ð Þ, g FeHS r ð Þ and g FeHM r ð Þ (SI, Figure S5).
At 21 GPa and 2,500 K, the number of hydrogen atoms in the metallic phase becomes larger than that in the silicate within 5 ps, but a few hydrogen atoms remain in the silicate over the entire simulation period of 20 ps (Figures 1; SI, Figure S4 and Movies S1-S4). This is also reflected in the development of the RDF (SI, Figure S5 for the MSH system), where the amplitude of the principal peak in g OHS r ð Þ decreases rapidly with time and the principal peak in g FeHS r ð Þ emerges concurrently as a result of H atoms moving from the silicate to iron. By contrast, the interaction between H M and Fe, and between H M and O does not change with time as H M atoms remain in the metallic phase. Nevertheless, the principal peak in g OHS r ð Þ remains noticeable even after 20 ps at 21 GPa and 2,500 K, suggesting a weakly siderophile character of hydrogen.
At the three other P-T conditions explored (40 GPa, 2,800 K; 75 GP, 3200 K and 140 GPa, 3,500 K; SI, Table S1), all hydrogen atoms in the silicate phase move into the metallic phase within 10 ps and stay there.
By also monitoring the distribution of hydrogen in the systems containing carbon (MSHC) and oxygen (MSHO) in the metallic liquid (Figures 1; SI, Movie S1), the two-phase results show that the siderophile nature of hydrogen established for MSH does not change. The initially small interaction between O and H M identified in the g OHM r ð Þ disappears completely; similarly, the principal peak in g OHS r ð Þ vanishes with time.

H 2 O Volume in Silicate and Metallic Melts
The molar volume of water in the metal V H M is predicted smaller than that in the silicate V H S in the DFT-MD simulations at all three P-T conditions considered further in thermodynamic integration (Figures 2a-2c) which constitutes an important aspect in the stronger affinity of hydrogen to metallic liquids relative to silicate melts. The volume difference Δ r V for reaction 2 decreases significantly at low P (5-25 GPa range) for the 2,500 K DFT-MD simulations, from −36% 5 GPa to −8% at 20 GPa for Δ r V/V H S . This is due to the fact that the silicate melt structure undergoes radical atomic arrangement from four-to six-fold Si-O coordination (SI, Figure S6a): While the nearest neighbor O-H distance represented by the principal peak position in g HO (r) (SI, Figure S6b) remains at virtually the same distance, the second neighbor O-H distances are shifted to considerably lower r with increasing P, reflected in the decrease of V H S with P. By contrast, the average interatomic distance between Fe and H in the metallic melts remains constant over the same P-range (SI, Figure S6d), similar to the results reported in a previous DFT-MD study (Posner & Steinle-Neumann, 2019). At higher P, differences are Δ r V/V H S = −16% (40 GPa, 4,000 K) and −12% (130 GPa, 4,000 K), and they contribute more significantly to G by a PΔ r V energy component due to large P: The PΔ r V term increases from −0.13 eV at 20 GPa and 2,500 K, to −0.54 eV at 40 GPa and 4,000 K, and −0.88 eV at 130 GPa and 4,000 K (Figures 2d-2f; SI, Table S2).

Thermodynamic Integration
Thermodynamic parameters corresponding to different hydrogen partitioning conditions are summarized in Tables S2-S5 in the SI. Figure S7 in the SI shows the variation of the radial distribution function g SiO (r) for Figure 1. Initial and final configurations for the two-phase simulations of silicate and metal liquid in contact. The initial configurations for our simulations were generated by performing independent DFT-MD simulations on single-phase H-bearing silicate and metal liquids at comparable P and identical T (SI, Figure S1). Metal and silicate configurations were then combined into a single cell containing a total of (a) 347 atoms for Mg 37 Si 37 O 111 H 6 -Fe 150 H 6 (MSH system), (b) 377 atoms for the carbon-bearing system Mg 37 Si 37 O 111 H 6 -Fe 150 C 30 H 6 (MSHC system) and (c) 353 atoms for the oxygen-bearing Mg 37 Si 37 O 114 H 6 -Fe 150 O 3 H 6 (MSHO system). The simulation duration is 20-25 ps (SI , Table S1). Animations for the development of the cells at all P-T conditions can be found in the supporting online material (SI, Movie S1). the (MgSiO 3 ) 15 H 2 O liquid and g FeFe (r) for the Fe 50 H 2 O liquid with changes in λ at 20 GPa and 2,500 K. As λ decreases, the first peak in g(r) decreases in amplitude, becomes broader, and tends to shift to a smaller distance, to finally disappears at λ < 0.0097, reflecting ideal gas behavior. Potential energy differences between the ideal gas and DFT liquid as a function of λ are shown in Figure S8 in the SI, with a smoothly varying integrand of eq. 5 at λ ≥ 0.0097, and a strong increase for smaller λ. In order to perform the integration, a variable transformation λ is done (Dorner et al., 2018), and eq. 5 becomes This transformation is chosen such that the transformed functional f [λ(x)]λ(x) k can be set to 0 at x = −1. Using Gauss-Lobatto quadrature, the integral becomes a sum where we chose k = 0.8 and a node number of n = 8, again following Dorner et al. (2018).
Calculated Δ r G and hydrogen metal-silicate equilibrium constant (expressed as log 10 K D , eq. 4) are shown in Figure 3 and given in Table S6 in the SI. The predicted K D value at 20 GPa (2,500 K) agrees with partition experiments of Okuchi (1997) at 7.5 GPa, and is in reasonable agreement with the experiments of Clesi et al. (2018) at P~20 GPa (SI ,  Table S6), although our results suggest a slightly stronger siderophile behavior. Differences in partitioning of hydrogen between the experiments by Okuchi (1997), on the one hand, and Clesi et al. (2018) and Malavergne et al. (2019), on the other hand, stem from the fact that Okuchi (1997) accounts for bubble volume in the recovered metal that  Table S2). (d-f) P · Δ r V contribution to the Gibbs energy change in reaction 2.  reflects exsolved hydrogen during quench due to its low solubility in solid Fe at ambient conditions (Antonov et al., 1998(Antonov et al., , 2019Fukai & Suzuki, 1986). Hirose et al. (2019) and Umemoto and Hirose (2020) similarly arrive at the conclusion that the experiments by Clesi et al. (2018) and Malavergne et al. (2019) underestimate the hydrogen content in the metal.
As we increase P to 40 GPa and T to 4,000 K, our results predict that the vast majority of hydrogen partitions into the metallic melt, leaving only 2% hydrogen in the coexisting silicate melt. Zhang and Yin (2012), by contrast, predicted hydrogen to be slightly lithophile at comparable P-T conditions. This discrepancy is likely due to the large interface-to-volume ratio of the two-phase model in Zhang and Yin (2012), which could induce major artifacts (Hong & Van De Walle, 2013). At 130 GPa and 4,000 K, K D is predicted to be exceedingly large such that the amount of hydrogen in silicate melt is negligible.
One should bear in mind that Δ r G of exchanging an H 2 O composition between silicate and metallic melts in reaction 2 as listed in Table S6 of the SI is very small compared to those typically associated with chemical reactions: for the reaction (CO) 50 + (O 2 ) 25 → (CO 2 ) 50 in the gaseous phase, for example, Δ r G = 130 eV (at 298 K, 1 bar) (Atkins et al., 2010). This raises the question of uncertainties in calculating Δ r G (and consequently K D ). The biggest source of uncertainty is the statistical fluctuation of U in DFT-MD simulations (for all λ in the thermodynamic integration). An uncertainty of 5 eV (~2%) in G for each of the components in reaction 2 translates to 10 eV for Δ r G, with severe consequences for values of K D : at 20 GPa and 2,500 K, log 10 K D would vary between −19 and 21. In order to decrease the error, two approaches would be useful: i. Larger simulation cells, with the goal to increase the number of exchanged H 2 O, would decrease statistical fluctuations (increased precision). ii. Simulations for F on a larger set of V-T points, with the goal to achieve higher accuracy by fitting physically motivated EoS models on this set, following prior work on silicate (de Koker & Stixrude, 2009) and metal liquids (Vlček et al., 2012).
However, both of these approaches require computational resources that are beyond current capabilities. Nevertheless, an empirical evaluation of our results by computing Δ r G of reaction 2 over shorter run durations of the DFT-MD trajectory shows negligible changes that do not influence log 10 K D values.

A Hydrogen-Rich Earth's Core
Using two independent computational approaches within the general framework of DFT-MD, two-phase simulations and thermodynamic integration, we predict that hydrogen strongly partitions into liquid metal over silicate, with P leading to a strong increase in the equilibrium constant (Figure 3). In a planetary context, this implies that for a geochemically reasonable hydrogen content of the primitive mantle (~1,000 ppm H 2 O) (e.g., Palme & O'Neill, 2013;Rubie et al., 2015) and assuming metal-silicate equilibration at P~40 GPa (Li & Agee, 1996;Siebert et al., 2012) for single-stage core formation, the hydrogen content in the core is close to 1 wt% H. This amount of hydrogen, together with some amounts of other light elements (e.g., Si, cf. Section 4.2), is consistent with the requirement to match ρ and velocity of the Earth's core (Umemoto & Hirose, 2020). It is further helpful to maintain the outer core in the molten state at relatively low T (Nomura et al., 2014;Sakamaki et al., 2009). However, it remains contentious whether Earth accreted wet or dry. In the latter case, as favored by the late veneer hypothesis (e.g., Albarède, 2009;Righter, 2003;Wang & Becker, 2013), no significant amount of water was incorporated into the core during differentiation.
Hydrous minerals such as δ-AlOOH, phase-H MgSiO 2 (OH) 2 , and their solid solutions have been shown to be stable over the entire range of lower-mantle P-T conditions (Duan et al., 2018;Ohira et al., 2014;Sano et al., 2008;Walter et al., 2015), and could therefore transport water into the lowermost mantle in a subducting slab. With the steep T-increase in the D " layer just above the CMB, these hydrous phases are likely to undergo dehydration, leading to hydrous silicate melts. Previous DFT-MD simulations (Bajgain et al., 2015;Mookherjee et al., 2008) suggest that such hydrous silicate melts are neutrally or negatively buoyant and therefore gravitationally stable at the base of the lower mantle, inhibiting efficient upward transport of volatiles through the mantle. In such a scenario, given the large value of K D at CMB conditions (Figure 3), the core would serve as an infinite sink for any hydrogen that comes in contact with it.

Mg/Si Fractionation Due to Core Formation
Beyond the preferred partitioning of H into the metallic melt, our two-phase simulations reveal a substantially different distribution of Si and Mg (SI, Figure S9) with a significant number of Si atoms (e.g., 8-10/37 for 41 GPa and 2,800 K in the MSH system) stably present in the metal for the last 5 ps of the atomic trajectories, while all Mg atoms remain in the silicate portion of the simulation cell. The absence of Mg from the metallic melts is not surprising, as our simulations were conducted at relatively low T that prohibit extensive mixing of silicate and metallic phases (Arveson et al., 2019;Huang et al., 2019). As Mg and Si have comparable diffusion coefficients (SI, Figure S10 and Table S7), the difference in behavior cannot be associated with their dynamics. Rather, their varying solution in Fe has to be associated with the fact that Mg has a weaker chemical affinity to metallic liquids than Si. Chemical reactions between metal and silicate investigated by DAC experiments (Badro et al., 2016;Takafuji et al., 2005) have shown that Mg solubility in liquid iron is extremely low for equilibration below 3,000 K, while up to 3 wt% Si dissolve into liquid iron, supporting the inference of our simulations.
This suggests sequestration of Si over Mg from the magma ocean during core formation, which in turn could explain the difference in the Mg/Si ratio between the chondritic value of~1.0 (Allègre et al., 1995) and that of the bulk silicate Earth (~1.25) (McDonough & Sun, 1995), supporting the hypothesis of Si incorporation in the core (e.g., Georg et al., 2007;Shahar et al., 2009) over Si enrichment in the lower mantle (e.g., Murakami et al., 2012).

Conclusions
We present results from two-phase and thermodynamic integration simulations based on density functional theory molecular dynamics that model hydrogen partitioning between liquid silicate and metal. The two-phase simulations show strong accumulation of hydrogen in the metallic portion of the simulation cell that we quantify in terms of an equilibrium constant K D by thermodynamic integration. Despite significant uncertainty, we find strong evidence for a moderately siderophile behavior of hydrogen at low pressure (20 GPa and 2,500 K), in quantitative agreement with the experiments by Okuchi (1997) and in qualitative agreement with the experiments of Clesi et al. (2018) at~20 GPa. We further describe a strong increase in K D with pressure, suggesting that hydrogen is strongly enriched in the core-forming metal during magma ocean segregation.
The two-phase simulations indicate non-negligible transfer of Si from the silicate to the metallic liquid, while Mg completely remains in the silicate. This qualitative evidence provides a plausible scenario for Si incorporation into Earth's core, potentially explaining the higher Mg/Si ratio of the bulk silicate Earth over chondritic values. Acknowledgments L.Y. has been funded by the Bayerisches Geoinstitut for his postdoctoral research. G. S. N. was funded by the Deutsche Forschungsgemeinschaft (German Science Foundation, DFG) grant STE1105/12-1 in the focus program "Building a habitable Earth" (SPP 1833). Computations were performed at the Leibniz Supercomputing Center of the Bavarian Academy of Sciences and the Humanities. The Vienna Ab Initio Simulation Package (VASP) is proprietary software for which licenses are issued at https://www.vasp.at/. In compliance with the journal's FAIR data availability requirements, the simulation files are described in computational methods and are archived in Figshare (10.6084/m9. figshare.12086184).