https://orcid.org/0009-0000-6600-9583
SUMMARY
Absolute palaeointensities are notoriously hard to obtain, because conventional thermal Thellier palaeointensity experiments often have low success rates for volcanic samples. The thermal treatments necessary for these experiments potentially induce (magnetic) alteration in the samples, preventing a reliable palaeointensity estimate. These heating steps can be avoided by pseudo-Thellier measurements, where samples are demagnetized and remagnetized with alternating fields. However, pseudo-Thellier experiments intrinsically produce relative palaeointensities. Over the past years, attempts were made to calibrate pseudo-Thellier results into absolute palaeointensities for lavas by mapping laboratory induced anhysteretic remanent magnetizations (ARMs) to the thermally acquired natural remanent magnetizations (NRMs). Naturally occurring volcanic rocks, however, are assemblages of minerals differing in grain size, shape and chemistry. These different minerals all have their own characteristic mapping between ARMs and thermal NRMs. Here, we show that it is possible to find these characteristic mappings by unmixing the NRM demagnetization and the ARM acquisition curves into end-members, with an iterative method of non-negative matrix factorization. In turn, this end-member modelling approach (EMMA) allows for the calculation of absolute palaeointensities from pseudo-Thellier measurements. We tested our EMMA using a noise-free numerical data set, yielding a perfect reconstruction of the palaeointensities. When adding noise up to levels beyond what is expected in natural samples, the end-member model still produces the known palaeointensities well. In addition, we made a synthetic data set with natural volcanic samples from different volcanic edifices that were given a magnetization by heating and cooling them in a controlled magnetic field in the lab. The applied fields ranged between 10 and 70 µT. The average absolute difference between the calculated palaeointensity and the known lab field is around 10 µT for the models with 2–4 end-members, while the palaeointensity of almost all flows can be retrieved within a deviation of ±20 µT. The deviations between the palaeointensities and the known lab fields are almost Gaussian distributed around the expected values. Although the two data sets in our study show that there is potential for using this end-member modelling technique for finding absolute palaeointensities from pseudo-Thellier data, these synthetic data sets cannot be directly related to natural samples. Therefore, it is necessary to compile a data set of known palaeointensities from different volcanic sites that recently cooled in a known magnetic field to find the universal end-members in future studies.
1. INTRODUCTION
The Earth's magnetic field is generated by convection in the liquid outer core of our planet. Since the 1970s, this magnetic field has been measured continuously around the world by satellites. To understand how the geomagnetic field behaves over time and possibly predict its future, it is important to have an accurate understanding of its more distant past. This would, for example, provide boundary conditions for models describing the behaviour of the geodynamo (e.g. Aubert et al. 2013; Sprain et al. 2019; Meduri et al. 2021), and will benefit models of (regional) palaeosecular variation that can be used for dating of archaeological artefacts and volcanic products (e.g. Korte et al. 2011; Nilsson et al. 2014; Pavón-Carrasco et al.2014).
The Earth's magnetic field is recorded by iron-bearing minerals in volcanic products, often basaltic lavas, when they cool below their Curie temperature. These volcanic products, therefore, may provide spot readings of the past state of the Earth's magnetic field for their moment of cooling and location on the planet. Full-vector information of the Earth's magnetic field from these lavas consist of a direction of the magnetic field, a palaeodirection and the strength of the field, a palaeointensity. While palaeodirections are usually relatively easy to obtain if the material has not undergone movement since its cooling, palaeointensity measurements are much more complicated (Tauxe & Yamazaki 2015). In principle, it is possible to estimate the intensity of an ancient magnetic field because the mechanism by which volcanic rocks acquire their primary magnetization is in theory linearly related to the ambient field for low fields such as the Earth's (Folgeraiter 1899; Koningsberger 1938). An important assumption of all palaeointensity techniques is therefore that the natural remanent magnetization (NRM) of the samples is a pristine thermoremanent magnetization (TRM), acquired at the time of cooling.
The most common method for obtaining absolute palaeointensities from samples with a TRM, is the Thellier–Thellier method (Thellier & Thellier 1959) and variations of their initial protocol (e.g. Coe 1967; Aitken et al. 1988; Tauxe & Staudigel 2004; Yu et al. 2004). The concept of a Thellier-style experiment is to stepwise replace the NRM by a laboratory induced TRM. Then, the assumed linear relationship between the NRM and the laboratory induced TRM, since they are both thermally induced, can be used to obtain the palaeointensity with the help of an Arai plot (Arai 1963; Nagata et al. 1963). In practice, however, Thellier-style measurements are hampered in many ways, often arising from the complex mineralogy of the iron-bearing minerals. For instance, chemical alteration during the experiment, complex magnetic behaviour of large iron-oxide minerals (i.e. ‘multidomain behaviour), and magnetic anisotropy of the sample may all prevent a successful determination of the palaeointensity (Tauxe & Yamazaki 2015). Some of these problems arise due to the subsequent heating steps that are necessary for a Thellier-style experiment. Avoiding heating the samples may prevent chemical alteration and to some extend multidomain behaviour to occur (e.g. Yu et al. 2002a). A Thellier-based experimental technique that does not use heating is the ‘pseudo-Thellier’ technique (Tauxe et al. 1995), which instead uses alternating fields (AF) to remove the NRM and impart laboratory magnetizations, in this case anhysteretic remanent magnetizations (ARMs), in the samples.
The pseudo-Thellier Technique was originally developed to help determining relative palaeointensity records from sediment cores (Tauxe et al. 1995), but the rationale of using AF steps instead of heating was later also applied to volcanic rocks. The main challenge with substituting the heating steps for AF steps in a Thellier-style experiment is that TRM demagnetization and remagnetization behaviour is not necessarily (linearly) proportional to AF demagnetization and remagnetization behaviour of the same minerals. Therefore, the slope of the linear fit in a pseudo-Thellier Arai diagram cannot directly be used to calculate the palaeointensity, and a conversion or ‘mapping’ of the AF demagnetization and remagnetization behaviour to the thermal demagnetization and remagnetization behaviour is necessary. Yu et al. (2002b) studied the analogy between ARM and TRM behaviour in volcanic rocks. They found that although the behaviour of ARMs and TRMs in volcanic rocks show some resemblance, it is difficult to map the behaviour one-to-one. Moreover, the mapping between ARMs and TRMs differs for iron oxides that differ in chemistry, size and shape. Since the magnetic carriers in volcanic rocks consists of assemblages of iron oxides with differing properties, the mapping between the TRM and ARM behaviour of the entire sample may very well be complex and nonlinear.
de Groot et al. (2013a) applied a strict rock-magnetic selection criterion to choose a group of samples with relatively similar ARM and TRM behaviour for a pseudo-Thellier study subjecting Hawaiian lavas. This allowed establishing a calibration equation by using a group of recent (1840–2010 AD) lavas that cooled in a known magnetic field for pseudo-Thellier experiments. This calibration relation was defined as a linear relation between the slopes of pseudo-Thellier measurements and their known palaeointensities. Surprisingly, this linear calibration relation has a non-zero y-axis intercept. This is problematic because this implies that if there is no magnetization present when cooling, the sample still has a palaeointensity of around 15 μT. Furthermore, this calibration relation is only applicable to samples that pass a strict selection criterion, which means that it is not applicable to all samples. Paterson et al. (2016) addressed the issue of the non-zero intercept of the y-axis, by using samples that were given a known TRM in the laboratory before the pseudo-Thellier experiments were done and anchoring the calibration relation to the origin. Theoretically, this calibration relation is more correct, but unfortunately it does not reproduce known palaeointensities from recent lavas better than the experimentally estimated calibration relation of de Groot et al. (2013a, 2015, 2016). Previous efforts to define a calibration relation for pseudo-Thellier experiments, that is, a mapping between ARM and TRM demagnetization and remagnetization behaviour, are hampered by the fact that they attempt to provide a single mapping between the TRM and ARM behaviour for bulk samples that may consist of assemblages of iron oxides that differ in magnetic properties. Because of these differences in magnetic properties, different minerals in a bulk sample may exhibit different relations between ARM and TRM behaviour, while they are measured simultaneously. This adds complexity to finding a proper mapping between ARMs and TRMs in a sample.
Here, we apply an end-member modelling approach (EMMA) to describe the magnetization of bulk samples as a combination of several empirical end-members with their own characteristic mapping between their ARM and TRM behaviour. After these end-members are defined, it is possible to unmix the measurements of bulk samples and obtain absolute palaeointensities from their pseudo-Thellier results. Such an EMMA has been successfully used in palaeomagnetic studies before; it was used to unravel and characterize curves of Isothermal remanent magnetization (IRMs, Kruiver et al. 2001) and ARMs (Robertson & France 1994; Egli 2003; Heslop & Dillon 2007). Here, we expand the unmixing of Heslop & Dillon (2007) to optimize for the two data sets in a pseudo-Thellier study (NRM demagnetization and ARM acquisition) simultaneously with the ultimate aim of retrieving the absolute palaeointensities from lavas. To illustrate the potential of the EMMA technique for pseudo-Thellier experiments on lavas we first apply it to a numerically created synthetic data set. Then, we apply it to a data set of basalts that were given a full TRM in the laboratory prior to measuring the pseudo-Thellier experiments.
2. METHODS
A palaeomagnetic sample used for palaeointensity studies typically consists of a mixture of different magnetic minerals that differ in size, shape and chemistry. This implies that these magnetic minerals also differ in magnetic behaviour. The rationale of the EMMA pseudo-Thellier technique is to identify common magnetic components, that is, end-members, in AF demagnetization and ARM acquisition measurements of volcanic samples. These end-members can be obtained by unmixing a large data set of AF measurements from volcanic samples with a known palaeointensity. Once these common end-members have been obtained they can be used to unmix measurements from volcanic samples, with an unknown palaeointensity, as a combination of these end-members with defined behaviour to obtain palaeointensity estimates.
The unmixing modelling technique used in this paper largely follows Heslop & Dillon (2007). They presented a modelling procedure that uses linear combinations of end-members to represent coercivity spectra. Importantly, the procedure of Heslop & Dillon (2007) unmixes and finds end-members for only one data set; we will expand on this and make the routine suitable to unmix two data sets simultaneously.
In general, EMMA starts with the measured (magnetization) data, for which the stepwise measurements of all the samples are contained in data matrix X. This matrix is formed with all measurements of each sample in a row and the measurements of the same demagnetization steps in each column. The rationale is to unmix data matrix |${\rm{X}}$| into the product of two matrices, |${\rm{S}}$| and A. Therein, matrix S contains the end-member curves, the common magnetic components, and matrix |${\rm{A}}$| the (non-negative) abundances, which contains the number of end-members present in each sample and their relative abundance (eq. 4). In practice, the measurements contain errors, and so error |$\epsilon $| must be added, where epsilon denotes the additive measurement noise.
To find the combination of end-members and abundances that represent the data best, a squared Euclidean distance cost function needs to be minimized to estimate the end-members and abundances, matrix |${\rm{S}}$| and |${\rm{A}}$|. The cost function gives a measure of the difference between data |${\rm{X}}$| and the modelled data |${\rm{AS}}$| (eq. 5):
This cost function can be minimized with respect to |${\rm{A}}$| and |${\rm{S}}$| (eqs 6 and 7), by using an iterative method of non-negative matrix factorization (Lee & Seung 2000). This is an optimalization problem, with the aim to find the best-fitting solution, that is, a combination of matrices A and S that explain the data X best.
Here, we aim to find the absolute palaeointensity from pseudo-Thellier measurements, which consists of two data sets; NRM demagnetization and ARM acquisition data. Therefore, we will unmix two data sets simultaneously by building on concept of unmixing only one data set introduced by Heslop & Dillon (2007) (eq. 4), namely:
where |${\rm{X}}$| and |${\rm{Y}}$| represent the NRM demagnetization and the ARM acquisition data sets, respectively. The absolute palaeointensity and laboratory field strength are indicated by |${{\rm{B}}}_{{\rm{ancient}}}$| and |${{\rm{B}}}_{{\rm{lab}}}$|. Both data sets have their own set of end-members, |${{\rm{S}}}_{\rm{X}}$| and |${{\rm{S}}}_{\rm{Y}}$|, which are different but correlated to each other, since it is the expression of the same magnetic component but in a different measurement. Lastly, both data sets have the same matrix of abundances |${\rm{A}}$|, since the abundance of the magnetic components of a sample is the same in the AF demagnetization as in the ARM acquisition measurement.
Expanding EMMA from one to two data-series requires changes to eqs (5)–(7). This is done in two steps: step 1 (calibration) aims to estimate the common end-members, |${{\rm{S}}}_{\rm{X}}$| and |${{\rm{S}}}_{\rm{Y}}$|, from a data set that acquired its magnetizations in a known magnetic field (in other words |${{\rm{B}}}_{{\rm{ancient}}}$| is known). In step 2 (classification), we use the end-members obtained in step 1, to check how well they are able to produce absolute palaeointensities (|${{\rm{B}}}_{{\rm{ancient}}}$|) from a data set that acquired its magnetizations in an unknown magnetic field.
2.1 Step 1: calibration
We measure |$X$| and |$Y$|, and |${B}_{\mathrm{ ancient}}$| and |${B}_{\mathrm{ lab}}$| are known |$ \to $| we estimate |${S}_X$| and |${S}_Y$| (to find them we also need to find |$A$|).
In a very similar way to eq. (5), we obtain two Euclidean distance cost functions from eqs (8) and (9), which we must minimize simultaneously with respect to |${\rm{A}}$|, |${{\rm{S}}}_{\rm{X}}$| and |${{\rm{S}}}_{\rm{Y}}$|. Note that in this equation |${\rm{X}}$| and |${\rm{Y}}$| have been divided by the known palaeointensities |${{\rm{B}}}_{{\rm{ancient}}}$| and |${{\rm{B}}}_{{\rm{lab}}}$|.
By using the same iterative method of non-negative matrix factorization as Lee and Seung (2001), the following equations for |${\rm{A}}$|, |${{\rm{S}}}_{\rm{X}}$| and |${{\rm{S}}}_{\rm{Y}}$| can be obtained:
where |${{\rm{\alpha }}}_1$| and |${{\rm{\alpha }}}_2$| are weighting factors which together sum to 1. These weighting factors determine the dependency on the AF demagnetization data set with respect to the ARM acquisition data set for the calculation of matrix |${\rm{A}}$|. These weighting factors can be empirically tested, to see which values give the most accurate calculation of the palaeointensity. The unmixing scheme starts with an initial guess for each end-member and the abundances, which should be slightly different for each end-member, otherwise the unmixing does not start properly. Any adequate initial guess gives similar results, observed from testing different initial guesses for the end-members.
2.2 Step 2: classification
We measure |$X$| and |$Y$|, and |${S}_X$|, |${S}_Y$| and |${B}_{lab}$| are known |$ \to $| we estimate |${B}_{ancien{\rm{t}}}$| (to find |${B}_{ancient}$| we also need to find |$A$|)
The aim of step 2 is to find the absolute palaeointensity (|${{\rm{B}}}_{{\rm{ancient}}}$|) for each pseudo-Thellier measurement from the common end-members that were defined in step 1. This can easily be obtained from eqs (8) and (9), since we have two equations and two unknowns, namely; |${\rm{A}}$| and |${{\rm{B}}}_{{\rm{ancient}}}$|.
2.3 Data sets
To illustrate the potential of EMMA for pseudo-Thellier experiments we use two different data sets in our study: a numerical data set and a synthetic data set based on basaltic samples that were given a full TRM in the laboratory. First, for 90 per cent of the data set (the training data set), step 1 is performed to find the end-members for that specific training data set. Second, step 2 is performed on the remaining 10 per cent of the data set (the test data set) to see how well absolute palaeointensities can be found from the end-members obtained in step 1. These two steps are repeated 100 times to get an average for the absolute calculated palaeointensities and an uncertainty estimate. For the EMMA technique, the data sets need to be normalized: all measurements, both the NRM demagnetization and ARM acquisition measurements, are normalized by the maximum value of the corresponding ARM measurement. In addition, during step 1, the data set is divided by the field strength of the ancient field (|${{\rm{B}}}_{{\rm{ancient}}}$|) in case of the NRM measurements and by the laboratory field strength (|${{\rm{B}}}_{{\rm{lab}}}$|) for the ARM acquisition measurement.
The first data set to be tested is numerical, to check whether EMMA works for optimizing two data-series at once. This numerical data set is constructed from a set of end-members (|${{\rm{S}}}_{\rm{X}}$| and |${{\rm{S}}}_{\rm{Y}}$|), chosen from the results of an arbitrary run of the model with the synthetic data set described below, which is multiplied by a randomly chosen |${\rm{A}}$| and |${{\rm{B}}}_{{\rm{ancient}}}$| to obtain the numerical data-series |${\rm{X}}$| and |${\rm{Y}}$|. In addition, we observe how well the model is able to handle this numerical data set with increasing amounts of Gaussian noise added to the data-series |${\rm{X}}$| and |${\rm{Y}}$|. This is done by defining the standard deviation of the Gaussian noise that is added to the data set by a percentage, |${\rm{\sigma \ error\ }}\mathrm{ per}\ \mathrm{ cent}$|, (2, 4, 6, 8 and 10 per cent) of the norm of the corresponding data set.
The second data set is synthetic and contains real pseudo-Thellier measurements of laboratory magnetized samples. The data set contains a total of 247 cores that were taken from five different volcanic edifices, namely the island of Réunion (France), Hawaii (USA), and Pico (Azores, Portugal); Mt. Etna on Sicily (Italy) and Iceland. The cores were remagnetized in an oven going up to 700 °C and were given a palaeointensity (|${{\rm{B}}}_{{\rm{ancient}}}$|) between 10 and 70 µT. The samples from each volcanic location were divided into six sets of about eight samples, while a set from each volcanic location was heated and cooled in one of the following field strengths (|${{\rm{B}}}_{{\rm{ancient}}}$|): 10, 22, 34, 46, 58 or 70 µT. Lastly, to check if the end-members of the synthetic data set have a physical meaning, that is, correspond to certain magnetic minerals and/or magnetic behaviour, a Curie temperature measurement was performed on a specimen from each set of samples.
3. RESULTS
To determine how well the EMMA technique is able to find the common end-members and calculate absolute palaeointensities for pseudo-Thellier measurements, first the numerical data set is used. The accuracy of the different models is described by two parameters. First, we use the difference in palaeointensity, |${{\rm{\Delta }}}_{{\rm{int}}}$|, which is defined as the average difference between the calculated palaeointensity by EMMA and the reference palaeointensity, for all samples; that is, |${{\rm{\Delta }}}_{{\rm{int}}}$| shows a potential systematic bias of the calculated palaeointensities from the model. Second, we define the absolute difference in palaeointensity |${| {\rm{\Delta }} |}_{{\rm{int}}}$| which is calculated in the same way as |${{\rm{\Delta }}}_{{\rm{int}}}$|, but takes the absolute values of the difference between calculated palaeointensity and reference palaeointensity, before calculating the average of all samples. This parameter shows on average how well the model is able to calculate the palaeointensities for different samples. When EMMA performs well, both these values should be close to zero. If the performance of EMMA breaks down, it is important to know whether the errors in predicted palaeointensities are random, as would be indicated by a near zero |${{\rm{\Delta }}}_{{\rm{int}}}$|, but a high |${| {\rm{\Delta }} |}_{{\rm{int}}}$|; or whether there is a bias in the results, in which case |${{\rm{\Delta }}}_{{\rm{int}}}$| would be negative for a systematic underestimate of the palaeointensity calculated by EMMA, or a positive |${{\rm{\Delta }}}_{{\rm{int}}}$| would show a systematic overestimate of the palaeointensity.
3.1 Numerical data set
For the numerical data set without added noise, EMMA is able to retrieve the palaeointensities perfectly: both |${{\rm{\Delta }}}_{{\rm{int}}}$| and |${| {\rm{\Delta }} |}_{{\rm{int}}}$| are close to zero (Fig 1). When adding Gaussian error to the numerical data set the calculation of the palaeointensity gradually becomes less accurate; |${| {\rm{\Delta }} |}_{{\rm{int}}}$| becomes linearly higher with increasing noise; |${{\rm{\Delta }}}_{{\rm{int}}}$|, which indicates systemic bias in the results is still very close to zero for an |${\rm{\sigma \ error\ }}per\ cent$| of 2 and 4, but then slowly decreases to negative values. It is important to note, however, that when adding error with a standard deviation of 10 per cent of the original data set, the model is still able to calculate the palaeointensities with an average absolute difference (|${| {\rm{\Delta }} |}_{{\rm{int}}}$|) of |$4.9{\rm{\ \mu T}}$|, which is better than the uncertainties in many palaeointensity studies on natural samples.
The performance of the 3 end-member model for a numerical data set with 255 measurements and |${{{\bf \alpha }}}_1 = \ {{{\bf \alpha }}}_2 = 0.5$|: illustrated by the average difference between calculated palaeointensity and reference palaeointensity for all samples (|${{{\bf \Delta }}}_{{{\bf int}}}$|) and the average difference in absolute calculated palaeointensity (|${| {{\bf \Delta }} |}_{{{\bf int}}}$|) plotted against an increasing amount of added error, averaged over 100 iterations. Where the standard deviation of the Gaussian noise added to the data set is determined by a percentage (2, 4, 6, 8 and 10 per cent) of the norm of the data matrix (X and Y, respectively).
3.2 Synthetic data set
The synthetic data set uses pseudo-Thellier results from 247 samples. The entire data set, consisting of an NRM demagnetization measurement and an ARM acquisition measurement for each sample, is normalized by the corresponding last ARM measurement. This normalization means, for example, that a maximum NRM value for a sample of 5, implies a five time as high NRM demagnetization measurement compared to its corresponding ARM acquisition measurement. Fig. 2 displays the distribution of maximum NRM values of the synthetic data set, normalized by the corresponding maximum ARM acquisition value. This shows that the bulk of the data set (89 per cent) has a maximum normalized NRM value between 0 and 9. The calculation of the palaeointensity becomes twice as accurate by removing the outliers which have a maximum normalized NRM value higher than 9. This choice is validated by the observation that measurements which have a maximum NRM value that is much higher than the maximum ARM values exhibit anomalous shapes. In addition, these extreme outliers are very dominant in an unmixing scheme, because EMMA is trying to find common end-members into which all the different NRM demagnetization measurements can be unmixed. For example, a very strong NRM demagnetization plot will require a strong end-member, which can be very disturbing for finding the common end-members which make up the bulk of the measurements. The difference in the NRM demagnetization data set (X) and ARM acquisition data set (Y) with and without the outliers is shown in Fig. 2.
(a)The distribution of the maximum NRM values, obtained from the first measurement of each NRM demagnetization measurement normalized by the last measurement of the corresponding ARM acquisition measurement. (b)a nd (c) The synthetic data set (247 samples), NRM demagnetization (X) and ARM acquisition (Y) data. (d) and (e) The data set (225 samples) selected from the distribution of the normalized maximum NRM values, the data with and normalized maximum NRM value between 0 and 9.
Both the numerical and synthetic data set give the most accurate palaeointensity calculation for an equal contribution of both the NRM demagnetization and ARM acquisition data sets for the calculation of A, in other words |${{\rm{\alpha }}}_1 = {\rm{\ }}{{\rm{\alpha }}}_2 = 0.5$|. For each iteration of the 3 end-member model, the average difference in palaeointensity (|${{\rm{\Delta }}}_{{\rm{int}}})$| and the absolute difference in palaeointensity (|${| {\rm{\Delta }} |}_{{\rm{int}}}$|) is calculated (Fig. 3a). In addition, the variance is calculated for each iteration to observe how well the calculated end-members (|${{\rm{S}}}_{\rm{X}}$| and |${{\rm{S}}}_{\rm{Y}}$|) and abundances (|${\rm{A}}$|) are able to explain the data matrices (|${\rm{X}}$| and |${\rm{Y}}$|), which EMMA is unmixing. The NRM demagnetization variance (|${\rm{NR}}{{\rm{M}}}_{\rm{v}}$|) and the ARM acquisition variance (|${\rm{AR}}{{\rm{M}}}_{\rm{v}}$|) are calculated by the difference between the original |${\rm{X}}$| and |${\rm{Y}}$| matrices and the |${\rm{X}}$| and |${\rm{Y}}$| matrices calculated from the unmixed end-members and abundances (Fig. 3b). It becomes clear from the values per iteration in Fig. 3 that the model quickly converges to a (local) minimum, with |${{\rm{\Delta }}}_{{\rm{int}}}$|, |${| {\rm{\Delta }} |}_{{\rm{int}}}$|, |${\rm{NR}}{{\rm{M}}}_{\rm{v}}$| and |${\rm{AR}}{{\rm{M}}}_{\rm{v}}$| all moving towards fairly low values.
(a) Difference in palaeointensity (|${{{\bf \Delta }}}_{{{\bf int}}}$|) and absolute difference in palaeointensity (|${| {{\bf \Delta }} |}_{{{\bf int}}}$|) for the 3 end-member model of the synthetic data set for the first 5000 iterations of step 1, the unmixing. (b) The NRM demagnetization (|${{\bf NR}}{{{\bf M}}}_{{\bf v}}$|) and ARM acquisition error (|${{\bf AR}}{{{\bf M}}}_{{\bf v}}$|) calculated, for the first 5000 iterations of the 3 end-member model, by the variance between the original |${{\bf X}}$| and |${{\bf Y}}$| matrices and the |${{\bf X}}$| and |${{\bf Y}}$| matrices calculated from the unmixed end-members (|${{{\bf S}}}_{{\bf X}}$| and |${{{\bf S}}}_{{\bf Y}}$|) and abundances (|${{\bf A}}$|).
The synthetic data set can be unmixed into different numbers of end-members. The unmixed end-members for the 2, 3 and 4 end-member model of the synthetic data set are in Fig. 4. The (absolute) average error for the subsequent calculation of the palaeointensity (|${{\rm{\Delta }}}_{{\rm{int}}}$| and |${| {\rm{\Delta }} |}_{{\rm{int}}}$|) for the test data set after 100 iterations, is in Table 1. It is evident that for every end-member model, |${{\rm{\Delta }}}_{{\rm{int}}}$| shows a slight bias towards negative values and |${| {\rm{\Delta }} |}_{{\rm{int}}}$| is around |$10{\rm{\ \mu T}}$|. Fig. 5, which displays the error distribution of |${| {\rm{\Delta }} |}_{{\rm{int}}}$| for each sample of the 3 end-member model, illustrates that the palaeointensity of almost all flows can be retrieved within a deviation of ± |$20{\rm{\ \mu T}}$|, where most flows can be found with a much smaller error.
The end-members, for the NRM demagnetization and ARM acquisition data set, calculated in step 1, unmixing, for the synthetic data set for the(a) and (b) 2, (c) and (d) 3, and (e) and (f) 4 end-member model.
The distribution of the difference in calculated palaeointensity (|${| {{\bf \Delta }} |}_{{{\bf int}}}$|) of each sample in the test data set, averaged over 100 iterations, for the 3 end-member model of the synthetic data set.
Difference in calculated palaeointensity (|${{{\bf \Delta }}}_{{{\bf int}}}$|) and absolute difference in calculated palaeointensity (|${| {{\bf \Delta }} |}_{{{\bf int}}}$|) calculated for the test data set averaged over 100 iterations with standard deviation (σ) for the 2, 3 and 4 end-member model of the synthetic data set in |${{\bf \mu T}}$|.
| |${{\rm{\Delta }}}_{{\rm{int}}}$| (|${\rm{\mu T}}$|) | σ | |${| {\rm{\Delta }} |}_{{\rm{int}}}$| (|${\rm{\mu T}}$|) | σ | |
|---|---|---|---|---|
| 2 end-member model | −1.65 | 3.11 | 9.44 | 1.74 |
| 3 end-member model | −2.38 | 3.15 | 9.41 | 2.08 |
| 4 end-member model | −3.95 | 2.74 | 10.00 | 1.73 |
| |${{\rm{\Delta }}}_{{\rm{int}}}$| (|${\rm{\mu T}}$|) | σ | |${| {\rm{\Delta }} |}_{{\rm{int}}}$| (|${\rm{\mu T}}$|) | σ | |
|---|---|---|---|---|
| 2 end-member model | −1.65 | 3.11 | 9.44 | 1.74 |
| 3 end-member model | −2.38 | 3.15 | 9.41 | 2.08 |
| 4 end-member model | −3.95 | 2.74 | 10.00 | 1.73 |
Difference in calculated palaeointensity (|${{{\bf \Delta }}}_{{{\bf int}}}$|) and absolute difference in calculated palaeointensity (|${| {{\bf \Delta }} |}_{{{\bf int}}}$|) calculated for the test data set averaged over 100 iterations with standard deviation (σ) for the 2, 3 and 4 end-member model of the synthetic data set in |${{\bf \mu T}}$|.
| |${{\rm{\Delta }}}_{{\rm{int}}}$| (|${\rm{\mu T}}$|) | σ | |${| {\rm{\Delta }} |}_{{\rm{int}}}$| (|${\rm{\mu T}}$|) | σ | |
|---|---|---|---|---|
| 2 end-member model | −1.65 | 3.11 | 9.44 | 1.74 |
| 3 end-member model | −2.38 | 3.15 | 9.41 | 2.08 |
| 4 end-member model | −3.95 | 2.74 | 10.00 | 1.73 |
| |${{\rm{\Delta }}}_{{\rm{int}}}$| (|${\rm{\mu T}}$|) | σ | |${| {\rm{\Delta }} |}_{{\rm{int}}}$| (|${\rm{\mu T}}$|) | σ | |
|---|---|---|---|---|
| 2 end-member model | −1.65 | 3.11 | 9.44 | 1.74 |
| 3 end-member model | −2.38 | 3.15 | 9.41 | 2.08 |
| 4 end-member model | −3.95 | 2.74 | 10.00 | 1.73 |
4. DISCUSSION
4.1 Performance of EMMA in pseudo-Thellier studies
The results using the numerical data set illustrate that the method used is mathematically sound. Without any noise, EMMA produces the palaeointensities of all samples without error. The model is also robust against Gaussian error, up to percentages of noise which are higher than expected in natural samples. However, it must be noted that the Gaussian error added to the numerical model may not be entirely representative of natural noise.
To test how well EMMA performs on samples that resemble natural samples, we made a synthetic data set using natural samples. This synthetic data set contains measurements of lab-magnetized samples from different volcanic edifices. The EMMA method is able to retrieve most palaeointensities (94 per cent) within a deviation of ±|$20{\rm{\ \mu T}}$| for a 3 end-member model of the synthetic data set. Most calculated palaeointensities are more precise, with 84 per cent being within a deviation of ±|$15{\rm{\ \mu T}}$|, 64 per cent within ±|$10{\rm{\ \mu T}}$| and 32 per cent within ±|$5{\rm{\ \mu T}}$|. In other words, the EMMA method is at least able to give an indication of most palaeointensities. Besides looking at the precision of the calculated palaeointensities, it is also important that the method is capable of calculating palaeointensities for a wide range of field strengths. Fig. 6 shows the distribution of the measured/calculated palaeointensity against the laboratory induced palaeointensity of the synthetic data set, for each sample. It is clear that, although far from perfect, EMMA is able to roughly estimate the palaeointensity for a wide range of field strengths. The average palaeointensities of the groups that were given palaeointensities 10, 22, 46 and 58 |${\rm{\mu T}}$| are very accurate, with 10.7, 21.3, 44.0 and 57.5 |${\rm{\mu T}}$|, respectively. The average of the 34 and 70 µT palaeointensity sets, however, deviate more with 27.1 and 61.3 µT, respectively. The deviation of the 34 and 70 µT palaeointensity sets is most likely because of the distribution of data in the synthetic data set. Due to removing the outliers in the beginning, more samples have been removed from these sets than from the other groups of samples. What is furthermore important to note is that the trendline in Fig. 6 goes almost through the origin, the y-intercept is only 0.89 µT. This implies that EMMA overcomes the problem of the non-zero y-axis intercept of the calibration relations in de Groot et al. (2013a, 2015, 2016), and does not need to be forced through the origin as in Paterson et al. (2016). The slope of the fit through all points, however, is slightly low (0.91), where we would expect 1.
![The blue dots give the distribution of the laboratory given palaeointensity [μT] against the measured/calculated palaeointensity [μT] of the test data set over 100 iterations by EMMA. The orange dots are the average measured/calculated palaeointensity [μT] for the entire palaeointensity set of 10, 22, 34, 46, 58 and 70 μT and the green line is the trendline through all the data points.](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/gji/235/3/10.1093_gji_ggad385/1/m_ggad385fig6.jpeg?Expires=1747927689&Signature=b74OKoaG5FD7b6isEZuULNu3a7lOTrFymgk1wowvFYOF58XpiQ2wH45-ivWFIDVEeXUOin4he76MHz4U3d9sfsrsXCb5SDEJdeMFNGJ9LvEVf0EHKsy0YbI2vmYqYR0-fOfswaSkXwBGeExw5LespSgXDY0v~nhMJGMmNBZfE4vRE1ac0ORRdyYZ-qLo4X5m6FB3vfGnTHhMbksMs-nnPCwWjv1li9eO8ExPtN3q3qi3F-ZJ6crNx7R3LCX0VX9QTRbqJhzsOLVRcfn-5oKL2Q3895A0ntThbZYSaovoHcEntOWjt8qG6nH28n3X4J04mHh~yfSWDEDGsYpwurz2zA__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
The blue dots give the distribution of the laboratory given palaeointensity [μT] against the measured/calculated palaeointensity [μT] of the test data set over 100 iterations by EMMA. The orange dots are the average measured/calculated palaeointensity [μT] for the entire palaeointensity set of 10, 22, 34, 46, 58 and 70 μT and the green line is the trendline through all the data points.
4.2 Physical representation of the end-members
To assess whether the end-members as produced by EMMA have any physical meaning, that is, represent groups of grains with similar magnetic behaviour, we measured the Curie temperature of each set of samples. Fig. 7 shows the correlation between the average abundances of each sample group (|${\rm{A}}$|) and the dominant Curie temperature (|${{\rm{T}}}_{\rm{C}}$|) of the 3 end-member model in two steps over 100 iterations. In step 1, where the training data set is unmixed into end-members, and in step 2, when the palaeointensity of the test data set is calculated from the end-members. The relationship between the Curie temperatures and abundances are notably similar for steps 1 and 2, for the 3 end-member model (Fig. 7). This means that the distribution between end-members calculated in step 1, is the same distribution of end-members found when determining the palaeointensity in step 2. A second observation, which was not seen for the 2 and 4 end-member models, is a difference in the relationship between the Curie temperature and the abundance of each end-member. End-member 3 has higher abundances for low Curie temperatures (±110–270 °C), whilst end-member 2 has higher abundances for higher Curie temperatures (±480–580 °C). End-member 1 shows no clear relationship between the Curie temperature and abundances for both steps 1 and 2. It must, however, be noted that the relationship between Curie temperatures and abundances is a trend and certainly not a rule which is illustrated by flows that clearly deviate from this trend in both plots.
The correlation between the most dominant Curie temperatures (|${{{\bf T}}}_{{\bf C}}$|) and the average abundances per sample group of each end-member (|${{\bf A}}$|) for step 1, (a) calibration and step 2, (b) classification. A linear trendline has been fitted for each end-member.
4.3 Limitations of our model
In an ideal scenario, the end-members that are found by EMMA would have a physical meaning, that is, represent a specific type of mineral with well-defined magnetic behaviour. Based on the ternary diagram of common iron oxides in lavas (Fig. 8), a 4 end-member model should select ilmenite, ulvöspinel, magnetite and haematite. In practice, however, the Curie temperatures of end-members ilmenite and ulvöspinel are below room temperature and are therefore not distinguishable when measuring palaeomagnetic samples. Moreover, the variation in Curie temperatures, and therefore the magnetic behaviour of iron oxides in the samples, is a continuous spectrum governed by Ti-content and oxidation state. Therefore, unmixing the data set into a certain number of end-members still is a simplification of the system of iron oxides, and EMMA has to define end-members that represent somewhat random points in the ternary diagram. Lastly, the magnetic behaviour of iron oxides is not only determined by magnetic chemistry but also by grain size and shape, which makes a physical representation of end-members even more complex. Nevertheless, the end-members that were defined are somewhat related to magnetic minerology, as illustrated by the trend between abundance of end-members and the Curie temperatures of the sample (Fig. 7).
Ternary diagram showing the ulvôspinel–magnetite solid solution series with increasing temperature and the ilmenite–haematite solid solution (after Readman et al. 1972).
Another limitation of the unmixing technique is the large range in values between the NRM curves, which makes it difficult to unmix the entire data set into a few end-members. For this reason, we choose to remove strong outliers, measurements of samples with a very high NRM curve, before running the EMMA routine. Despite knowing that the ratio between the maximum NRM and ARM values is determined by the chosen laboratory field strength during the ARM measurements, which in our case is always |$40{\rm{\ \mu T}}$|. This inherently favours a certain range of palaeointensities for the data sets presented in this paper. When further developing EMMA, it would be interesting to use measurements with a larger range of laboratory field strengths used during the ARM measurement.
Lastly, the method currently has no parameter or check that validates the accuracy of the palaeointensity calculation. The method is able to give an indication of the palaeointensity for almost all samples in the synthetic data set. However, the range of error for these calculations in palaeointensity is large. For the future development of this method, it would be very useful if such a validity parameter can be added to the EMMA method.
5. CONCLUSION AND OUTLOOK
We have shown, despite some limitations, that by unmixing a training data set into end-members and using these end-members to calculate the palaeointensity for a test data set we are able to at least provide a reasonable, first-order, approximation of the reference palaeointensity (Fig. 5), at least for our numerical and synthetic data sets. Our approximation of palaeointensities using EMMA is also better than previous attempts to calibrate pseudo-Thellier results from lavas into absolute estimates of the palaeointensity (de Groot et al. 2013a, 2015, 2016; Paterson et al. 2016). It is also important to emphasize that the EMMA method does not exclude samples based on rock-magnetic properties before the analysis and is applicable to all volcanic samples measured.
Interpreting pseudo-Thellier data from lavas using EMMA improves both the accuracy and the amount of data that can be interpreted. The numerical and synthetic data set are the first steps into defining a set of end-members that can be generally applied to volcanic or basaltic rocks from different volcanic edifices. Nevertheless, we know that palaeointensity experiments on samples that were just given a full TRM in the laboratory prior to the palaeointensity experiments always perform better than natural samples (e.g. de Groot et al. 2013b). Therefore, the data sets presented here cannot be directly related to natural samples, although they show that there is potential for using this end-member modelling technique for finding an approximation of the absolute palaeointensities from pseudo-Thellier data. To define generally applicable end-members, it is necessary to compile a sample set of volcanics that recently cooled in known palaeointensities from different volcanic sites, that is, that cooled in the Earth's magnetic field with different field strengths.
ACKNOWLEDGEMENTS
We thank our colleagues at palaeomagnetic laboratory Fort Hoofddijk for their help, lively discussions and comradeship. This project is funded by the Dutch Science Foundation (NWO) VIDI grant VI.Vidi.192.047 to LVdG. LVdG conceived of the presented idea. LvG developed the theory, carried out the experiments, performed the computations and took the lead in writing the manuscript. LVdG and TvL provided critical feedback and helped shape the research, analysis and manuscript. LVdG supervised the project.
CODE AND DATA AVAILABILITY
The EMMA.py code and the laboratory created test data set can be found on https://doi.org/10.5281/zenodo.8367125.
