Paleomagnetic field reconstruction from mixtures of titanomagnetites

a Department of Earth Science and Engineering, Imperial College London, SW7 2AZ, UK b Instituto Dom Luiz, Faculdade de Ciências, Universidade de Lisboa, 1749-016 Lisboa, Portugal c School of Earth Sciences, University of Bristol, Wills Memorial Building, Queen’s Road, Bristol BS8 1RJ, UK d Lamont–Doherty Earth Observatory, Columbia University, Comer Geochemistry Building, PO Box 1000, Palisades, NY 10964-8000, USA

Stepwise thermal demagnetization and alternating field (AF) demagnetization are commonly used in paleomagnetic studies to isolate remanent magnetic components of different origins. The magnetically hardest, i.e. highest unblocking temperature/peak field component is often interpreted as the primary magnetization and magnetically softer components as subsequent remagnetizations due to geological events posterior to the formation of the rock, such as reheating or formation of new magnetic minerals. The correct interpretation of the sequence of the geological events such as tectonic rotations from paleomagnetic data often relies on correctly attributing the observed magnetic directions to the remanence carriers and acquisition mechanisms. Using a numerical model to simulate remanence acquisition and stepwise thermal and AF demagnetization experiments, we show that the presence of mixtures of different magnetic minerals, such as magnetite and titanomagnetites of varying titaniumcontent can have very significant effects on Zijderveld plots. In thermal demagnetization experiments a spurious third component at intermediate temperatures or a continuous curvature may arise from an overlap of the primary remanence with a subsequent thermal or viscous remagnetization carried by small-grained iron-rich magnetite and large-grained titanium-rich titanomagnetite. AF demagnetization plots of magnetic mixtures are even more complex: primary and secondary remanences carried by different minerals may appear as either three or four components in Zijderveld plots. During alternating field demagnetization the highest coercivity component is not necessarily equivalent to the primary remanence and does not necessarily correspond to the highest temperature component in an analogous thermal demagnetization experiment, i.e., the primary remanence direction cannot be recovered. The effects are shown to be due to the different responsiveness of magnetite and titanomagnetites towards viscous or thermoviscous remanence acquisition: remanent magnetizations with long acquisition times are more effectively recorded by titanium-poor minerals, while short acquisition times are equally well recorded by titanium-rich minerals. In demagnetization experiments on laboratory timescales, the relative contribution of two minerals to Zijderveld plots differs to the relative contribution of remanence acquisition over geological timescales, leading to overlapping components in Zijderveld plots. The model was also used to simulate paleointensity (ancient magnetic field intensity) experiments and it was found that the grain distribution affects the slope of Arai plots, but is negligible compared to the effect of the cooling rate of NRM acquisition. The simulations suggest that for slowly cooled rocks a cooling rate correction of up to 1.5 to 1.6 may be required depending on the mineralogy.

Introduction
Paleomagnetic observations continue to provide constraints on some of the most fundamental theories of the deep Earth structure, the dynamics of near surface processes and the evolution and * Corresponding author.
E-mail address: t.berndt13@imperial.ac.uk (T. Berndt). development of the geodynamo (Tarduno et al., 2015;Biggin et al., 2015;O'Rourke and Stevenson, 2016). Reliable interpretation of paleomagnetic data can only be achieved through correct identification of the natural remanent magnetization (NRM) components and their directions; we are usually, but not always, interested in the primary remanent magnetization's intensity and its direction carried by the magnetic minerals within rocks. Among the most common magnetic minerals occurring in rocks are both stoi- where titanium atoms substitute the iron atoms at varying proportions x (Dunlop and Özdemir, 1997). In nature, rocks do not always contain only a single type of magnetic mineral but may contain mixtures, for example of titanomagnetites of varying compositions. The grain-sizes of the magnetic mineral have been found to correlate with the titanium content in oceanic basalts (Zhou et al., 1997(Zhou et al., , 2000 and the process of exsolution can move titanium cations in the crystal lattice of the Fe 3−x Ti x O 4 , accumulating them in some places and depleting them in others, thereby effectively creating an amalgam of high titanium content titanomagnetite grains and pure magnetite or low-titanium content titanomagnetite grains (Dunlop and Özdemir, 1997). To correctly interpret paleomagnetic signals of natural rocks, it is important to understand the effect of such magnetic mineral mixtures on the paleomagnetic recording fidelity. We developed a numerical model to predict the behavior of titanomagnetite mixtures with respect to three of the most fundamental paleomagnetic studies: (1) directional analysis in thermal demagnetization experiments, (2) directional analysis in alternating field (AF) demagnetization experiments, and (3) Thellier-type paleointensity estimates (Thellier and Thellier, 1959).

Model
A numerical model has been built, that simulates an assembly f (x, V ) of titanomagnetites of different titanium content x and different grain volumes V . The model is built on Néel (1949) theory of single-domain (SD) magnetic particles. The evolution of normalized magnetic moment n (magnetic moment divided by the spontaneous magnetization) with time is given by the differential equation (Néel, 1949) where τ is the relaxation time and n eq is the value of the normalized magnetic moment in thermodynamic equilibrium. The relaxation time is given by where τ 0 is the atomic attempt time, which was set to be 10 −10 s in the model (Berndt et al., 2015), μ 0 is the vacuum permeability, k is the Boltzmann constant and H 0 is the applied magnetic field. The equilibrium magnetic moment is given by a Maxwell-Boltzmann distribution The spontaneous magnetization at high temperature is modeled using the analytical approximation (Dunlop and Özdemir, 1997) and the microscopic coercivity H K is calculated assuming that shape anisotropy dominates, for which using a shape anisotropy factor N. For titanomagnetites shape anisotropy and magnetocrystalline anisotropy are relatively weaker than for magnetite, but magnetostriction increases (Dunlop and Özdemir, 2007). For simplicity, however, we assume strongly elongated grains with dominant shape anisotropy for all titanium contents with a common value of N = 0.5 for all grains and alignment of their elongation axis with the field.
The titanium content is assumed to have two effects: (1) it lowers the Curie temperature T C , and (2) it reduces the roomtemperature spontaneous magnetization M s0 . The Curie temperature is modeled by the quadratic equation where the coefficients a = 280 and b = 500 were found from a least-squares fit to the data published by Dunlop and Özdemir (1997), and T C ,T M0 = 580 • C is the Curie temperature of magnetite.

VRM and TRM acquisition
The grain distribution is discretized by a matrix of 1000 volumes V between 10 −24 and 10 −21 m 3 (being equal to cubes of 10 to 100 nm), separated on a logarithmic scale, and 100 equally spaced Curie temperatures T C between 0 • C and 580 • C (corresponding to various different titanium compositions x according to eq. (7), for clarity we quote T C rather than x values in the diagrams). The magnetization of each of these grains can take on any magnetization value representing a large number of grains, and not just ±1, as for a single SD grain.
For viscous remanent magnetization (VRM) acquisition at a temperature T A , the equilibrium magnetizations n eq (eq. (4)) and the relaxation times (eq. (2) and (3)) are calculated for each grain set (V , T C ) and the resulting new magnetization state n ne w is calculated from eq. (1). Thermoremanent magnetization (TRM) acquisitions are simulated by repeatedly following this procedure for 2000 temperature steps T i , decreasing by small temperature steps T until room temperature is reached. Various scenarios of different combinations of acquired VRMs and TRMs at different times and temperatures were run. Generally, linear cooling was used, but for one case Newtonian cooling was used for a paleointensity scenario, as cooling rates are known to have a significant effect on paleointensities (Dodson and McClelland-Brown, 1980;Halgedahl et al., 1980).

Thermal demagnetization
Step-wise thermal demagnetization was simulated by repeatedly applying VRMs at successively higher temperatures in zero field. This simulates the time at which the sample is kept at a high temperature in a thermal demagnetizer. After each step, the total remanent magnetization vector is calculated, which is the sum the magnetization vectors n (V , T C ) of all different grain sets, and the total spontaneous magnetization is calculated by summing the product of M s , volume V and the grain distribution f (V , T C ).

AF demagnetization
AF demagnetization is modeled based on the simplified assumption that all grains with a coercivity H C less than the maximum amplitude H of the alternating field get demagnetized. The coercivity is given by where H K is given by eq. (6) and H q is the thermal fluctuation field given by (Néel, 1949) Using eq. (6), and approximating the time t as half the inverse of the frequency f of the AF field (Worm, 1998), this simplifies to The amplitude H is successively increased and at each step the remaining total magnetization is calculated by summing the magnetization vectors n (V , T C ).

Paleointensity
A series of Thellier-type paleointensity experiments (Thellier and Thellier, 1959) were simulated following the methodology of Coe (1967). First a TRM acquisition was simulated using either linear or Newtonian cooling. Then, Arai plots (Nagata et al., 1963) were produced by simulating demagnetization steps to temperatures T i by calculating the viscous decay in zero field during heating to T i at 1 K/s, holding the temperature for 10 min and cooling back to room temperature at 1 K/s, and calculating the remaining NRM. Each step was followed by the simulation of a heating in zero field at 1 K/s, followed by a VRM acquisition in a 30 μT field for 10 min at T i , representing the hold time in the furnace in field in a Thellier-type experiment, followed by a TRM acquisition on cooling from T i to room temperature at 1 K/s, representing the in-field-cooling of a Thellier-type experiment.

Scenarios and grain distributions
The model was used to simulate the magnetic behavior of a number of different grain distributions for a number of different remanence acquisition scenarios.

Remanence acquisition scenarios
To investigate the effect of mixtures of titanomagnetites on vector demagnetization plots, the following three extreme scenarios of remanence acquisition were used: 1. A primary full TRM acquisition over 1 h was simulated, followed by a perpendicular pTRM acquired at 100 • C over 1 h. This acquisition time represents fast-cooling submarine lavas (Bowles et al., 2005). 2. A primary full TRM acquisition over 100 ka was simulated, followed by a perpendicular pTRM acquired at 100 • C over 100 ka. This timescale is typical of slowly cooling intrusive rocks . 3. A primary full TRM acquisition over 100 ka is simulated, followed by a perpendicular VRM acquired at room temperature (20 • C) over 100 ka.
The ambient magnetic field was set to be H 0 = 30 μT. For the stepwise thermal demagnetization experiments, a hold time of 10 min was used. It was expected from theory that for pure magnetite, all of these scenarios yield two perpendicular magnetic components in thermal demagnetization Zijderveld (1967) plots. They are expected to only differ in the unblocking temperatures of the remagnetizations: According to Pullaiah et al. (1975), scenario 1 should have an unblocking temperature close to 100 • C, as the timescale of the thermal demagnetization is similar to the timescale of acquisition, scenario 2 should have an unblocking temperature of 206 • C (linear cooling over 100 ka is equivalent to 5300 yr at constant temperature according to York 1978aYork , 1978b and scenario 3 should unblock at 135 • C. For the AF demagnetization experiments, a frequency of f = 50 Hz was used. Throughout the following treatment, the primary full TRM will be referred to as the characteristic remanent magnetization (ChRM), while the secondary remagnetization will be referred to as either partial TRM (pTRM) or VRM, respectively. For the paleointensity experiments, four scenarios were simulated: 1. A full TRM acquired by linear cooling over 100 ka in a 30 μT field. 2. A full TRM acquired by linear cooling over 1 h in a 30 μT field. 3. A full TRM acquired by linear cooling over 1 h in a 30 μT field, followed by a VRM in the same direction acquired over 100 ka in the same field. This scenario is meant to test if viscous overprints of rocks formed during the Brunhes chron have an effect on paleointensity determinations. 4. A full TRM acquired by Newtonian cooling over 100 ka in a 30 μT field. In order to avoid the cooling process to take infinitely long, the ambient temperature was set to 15 • C and the cooling rate was chosen such that a target temperature of 20 • C is reached after 100 ka.

Grain distributions
The model has been run with a series of different grain size and composition distributions. The first distribution investigated is a bimodal grain distribution with a magnetite (T C = 580 • C) peak at (10 nm) 3 , and a secondary peak around (30 nm) 3 large titanomagnetite TM60 with a Curie temperature of 200 • C (Fig. 1a). Note the smallest pure magnetite grains are superparamagnetic at room temperature and the remanence is due to only the grains larger than the peak value of 10 nm. The TM60 amounts to approximately 22% of the volume of the magnetic material, but only to 8% of the magnetic intensity due to the lower M s0 of TM60 and only 2.3% of the total number of grains due to the larger grain size. The Curie temperature of TM60 can be seen in the simulated M s (T ) curve in Fig. 2, with a value of 220 • C determined by the maximum second derivative method (Ade-Hall et al., 1965).
The second distribution was a continuous grain distribution with pure magnetite and titanomagnetite with a Curie temperature close to room temperature as its end-members, with the mean grain size increasing with titanium content (Fig. 1b). This distribution was chosen, because distributions with larger titanium-rich and smaller titanium-poor titanomagnetite grains have been observed in nature (Zhou et al., 1997(Zhou et al., , 2000. Such distributions lead to M s (T ) curves that decrease steeply at low temperatures before slowly leveling off at high temperatures before vanishing at the Curie temperature of magnetite (Fig. 2). Simple methods such as the maximum second derivative method (Ade-Hall et al., 1965) do not allow to obtain much insight into the mineralogy in this case: doing so would yield an intermediate value between the Curie temperatures of the end-members magnetite and low-titanium titanomagnetite and would miss the fact that the M s curve is due to a mixture of various minerals with a wide range of Curie temperatures.
The final case was a broad grain distribution that includes grains of all sizes and all titanium contents, with a slight correlation between grain-size and titanium-content, and log-normally distributed grain-volumes (Fig. 1c). The spontaneous magnetization curve shows a similarly sharp decay as the continuous magnetite-  titanomagnetite distribution at low temperatures, leveling off at high temperatures.
For the Thellier-type paleointensity experiments, additionally a pure magnetite grain distribution was investigated that equaled the distribution of the grains with T C = 580 • C of the bimodal grain distribution above.

Grain distribution 1: bimodal distribution
The stepwise thermal demagnetization plot of scenario 1 (Fig. 3a), reconstructs the directions of the two magnetic components as expected. Scenario 2, however, shows a demagnetization plot that could incorrectly be interpreted as having three magnetic components: one unblocking around 130-140 • C, with the expected direction of the pTRM, one intermediate direction unblocking around 200-210 • C and the original ChRM. The middle component is an artifact of two different magnetic minerals, however: A pTRM acquired at 100 • C during cooling over 100 ka should be removed at 206 • C in a 10 min demagnetization experiment, whereas the same pTRM acquired by titanomagnetite with a Curie temperature of 200 • C should be demagnetized at 133 • C (York, 1978a(York, , 1978bDodson and McClelland-Brown, 1980): The first (lowest temperature) apparent direction below 133 • C in the demagnetization plot corresponds to the demagnetization of the pTRM carried by the titanomagnetite, the second apparent direction up to 206 • C corresponds to the simultaneous demagnetization of the pTRM carried by magnetite and the ChRM carried by the titanomagnetite, and the third apparent direction above 206 • C corresponds to the demagnetization the ChRM carried (mostly) by the magnetite. A similar effect occurs in scenario 3 (Fig. 3g) with a VRM acquired over 100 ka at room temperature: up to 80 • C the direction of the pTRM is observed, above that a curvature up to 140 • C is seen and at higher temperatures the ChRM is recovered. In a 10 min demagnetization experiment, the unblocking temperature of such a VRM is 139 • C for magnetite and 79 • C for the titanomagnetite (Pullaiah et al., 1975). The curvature between these two temperatures is due to the overlap of the VRM carried by magnetite and the ChRM carried by titanomagnetite.
The AF demagnetization plots (Fig. 4a, 4d and 4g) all show four apparent components: the direction of the overprint is seen at low field H , after that a curvature approaching the ChRM direction is visible, followed by another section of the overprint's direction, and at highest fields the ChRM direction is again observed.

Grain distribution 2: continuous distribution
As in the case of a bimodal distribution, both magnetic components are accurately recovered for scenario 1 (Fig. 3b) in stepwise thermal demagnetization. Scenario 2 has three apparent directions, the first of which unblocks around 115 • C, and the second of which unblocks around 180 • C (Fig. 3e). These two unblocking temperatures of the pTRM correspond to titanomagnetite with Curie temperatures of 140 • C and 360 • C, respectively, for a demagnetization time of 10 min. Neither of these are the end-members of the titanomagnetite grain distribution (Fig. 2). While the first inflection point at 115 • C is relatively clear, the second one around 180 • C is curved, slowly approaching the final ChRM direction. Scenario 3 shows a similar trend with a lower blocking temperature of 45 • C and a higher one of 110 • C, which for the VRM of 100 ka, corresponds to titanomagnetite of Curie temperatures of 80 • C and 370 • C (Fig. 3h). A notable difference to the bimodal grain distribution is that the intensity of the VRM is significantly weaker, although the intensity of the pTRM is of a similar order. Hence this grain distribution is less responsive to VRM acquisition than the bimodal distribution (Fig. 3g).
The AF demagnetization plot of scenarios 1 and 3 ( Fig. 4b  and 4h) show similar trends to the ones of the bimodal grain distribution (Fig. 4a and 4g), but with different intensities and demagnetizing H fields: at low fields the pTRM/VRM direction is seen, but for higher fields an "S"-shape is observed, starting at approximately the ChRM direction, then bending into an intermediate direction and then bending back into the ChRM direction. The ChRM-direction can only be approximately isolated at the highest fields (>55 mT): the curvature due to overlap with the pTRM/VRM is small at the highest fields, it does not, however, necessarily vanish, such that the obtained ChRM direction may be imperfect. Compared to the bimodal distribution, the Sshape is greatly reduced in intensity in the continuous distribution and is dominated by the ChRM component. Additionally, like in the thermal demagnetization case (Fig. 3b and 3h) the VRM is strongly suppressed compared to the pTRM, almost disappearing completely at low coercivities. Therefore, in this particular grain distribution ( Fig. 4b and 4h), the ChRM-direction would be better isolated at lower fields (20-30 mT), contrary to the bimodal distribution ( Fig. 4a and 4g), where low fields (30-40 mT) showed an intermediate direction between the ChRM and pTRM/VRM directions. Scenario 2 (Fig. 4e) differs from the others in that the ChRM cannot recovered: after the VRM is isolated at coercivities below 23 mT, an S-shaped curve begins that approximates the ChRM direction between ∼30-40 mT, but turns into an intermediate direction at higher coercivities. No part of the diagram completely isolates the ChRM direction. Moreover, at the highest coercivities, the observed direction differs more from the ChRM direction than in the 30-40 mT range: the common assumption that the ChRM is best isolated at the highest demagnetization step is invalid.

Grain distribution 3: broad grain distribution
The demagnetization plots, both stepwise thermal demagnetization and AF demagnetization, all show a strong curvature between the ChRM and the pTRM/VRM (except scenario 1 in thermal demagnetization, Fig. 3c), contrary to the previous cases, where more than two distinct apparent components were observed. The curvature appears in a similar temperature range as in the case of the continuous magnetite-titanomagnetite distribution (Fig. 4), with the curvature lying between 110 • C and 190 • C for the pTRM (scenario 2) and between 50 • C and 120 • C for the VRM (scenario 3). Again, the VRM appears slightly weaker than the pTRM due to the presence of titanomagnetites that are less responsive to VRM acquisition.
The AF demagnetization plots do not show three to four apparent components as in the previous cases ( Fig. 3 and 4), but rather show a strong curvature between the pTRM/VRM and the ChRM. For this grain distribution, both the pTRM/VRM and the ChRM directions can be recovered from both stepwise thermal and AF demagnetization plots.

Paleointensity experiments
Arai plots (Nagata et al., 1963) were calculated for the three grain distributions and for a pure magnetite using the four scenarios described in section 3 (Fig. 5). All grain distributions show two types of behavior: a slow cooling behavior for linear and Newtonian cooling over 100 ka, and a fast cooling behavior for linear cooling over 1 h (with or without subsequent VRM acquisition). Within these two categories, Arai plots are almost identical for all samples, with the only exception of the 1 h TRM followed by a 100 ka VRM for the pure magnetite simulation (and, to a lesser degree, the bimodal distribution). This observation suggests that for paleointensity experiments, the effect of magnetic mineralogical mixtures is almost negligible, and the dominant factor impacting the slope of Arai plots is the cooling rate.
Arai plots are almost linear over the whole temperature range, but slight variations of the slope dM pT R M /dM N RM occur (Fig. 6).
The 1 h cooling scenario (circles) shows the most constant slope; slightly more than unity. A slope of one is expected if the NRM acquisition time equals the pTRM acquisition times in the Thellier/Coe-type experiment, but an exact comparison of the two timescales is difficult, as NRM acquisition occurs during cooling, while pTRM acquisition occurs during a hold time (i.e. VRM acquisition) of 10 min at elevated temperature and the subsequent cooling. The scenario also shows that the slope tends to increase slightly at higher temperatures depending on the grain distribution.
In contrast, the slow cooling scenarios show a much steeper slope around 1.4-1.6 (Fig. 6). Newtonian cooling (triangles) tends to further increase slopes at intermediate temperatures compared to linear cooling (squares), but only marginally. An obvious feature of the slow cooling scenarios is the strong increase of slope at low temperatures before reaching a peak around 100 • C and then slowly decreasing.
The scenario of a fast cooling TRM and a subsequent 100 ka VRM combines features of both the fast and slow cooling scenarios: at low temperatures, before the VRM is unblocked, the slope equals that of the slow cooling 100 ka TRM scenarios. At higher temperatures, the slope quickly approaches that of the fast cooling 1 h TRM scenario. This behavior is expected as the two parts of the Arai plots show two distinct magnetic components, that only coincide in direction, but not in intensity.

Thermal demagnetization
With the exception of the thermal demagnetization of the first scenario (fast-cooling TRM acquisition), all the simulated Zijderveld plots significantly deviate from the expected two-component behavior. This can be explained by considering the individual sets of grains that carry the remanence and who they get demagnetized.  (Nagata et al., 1963) in Fig. 5 as a function of temperature for different distributions and different remanence acquisition scenarios. Plots are normalized by maximum pTRM (normalization by pTRM rather than NRM was done because the pTRM was independent of the remanence acquisition scenario). First data point has non-zero pTRM as it corresponds to a "heating" to room-temperature, i.e., a VRM acquisition over 10 min.
The grains carrying the two remanent magnetizations are indicated in Fig. 7a, 7c and 7e for the bimodal distribution. In all three scenarios, the magnetite distribution as well as the titanomagnetite both partially carry the ChRM and partially the VRM/pTRM. The line separating the two magnetic components (dashed line) depends on both the acquisition temperature and the acquisition time; the effective acquisition time in the case of the pTRMs. All grains to the bottom left of the dashed line carry the remagnetization, whereas all grains to the top right of the line preserve the ChRM. Such a line can be calculated for any time and temperature. In general, increasing either the temperature or the time shifts the line to the top right. When demagnetizing the sample, the solid line is swept from the bottom left corner of the diagram to the top right corner as the temperature is increased, demagnetizing the grains below. The remaining remanence is carried by the grains to the top right of the line and measured after each heating. If both the time and the temperature of the demagnetization experiment are identical to the acquisition (as in scenario 1, Fig. 7a), both the acquisition and the demagnetizing lines are identical, but if the timescale of the demagnetization experiment differs to the acquisition timescale, then the slope of the demagnetizing line will differ to the slope of the acquisition line (Fig. 7c). While either one, increasing temperature or increasing time shifts the line to the top right, increasing time does so while tilting it clockwise, whereas increasing temperature does so while tilting it anti-clockwise. For this reason, the titanomagnetite is more responsive to increases in temperature, while the magnetite is more responsive to increases in time. As thermal demagnetization is usually done on a shorter timescale than acquisition (minutes to hours versus days to thousands of years), the titanomagnetite tends to be demagnetized first. As the demagnetization progresses, the larger titanomagnetite grains that preserve the ChRM and the smaller magnetite grains that carry the VRM/pTRM are demagnetized simultaneously and their components overlap. On further heating the larger magnetite grains that carry the ChRM are demagnetized. A similar effect occurs in the continuous grain distribution (Fig. 7b, 7d and 7f). Due to the strong correlation of grain volume with titanium content, the remagnetization (pTRM and VRM) affects two distinct grain populations: small-grained low-titanium magnetite as well as large-grained high-titanium content grains. In the first scenario, the same populations are activated during demagnetization as during acquisition, which is expected for any grain distribution if the acquisition time equals the demagnetization time. Therefore, the demagnetization plot shows only two components (Fig. 3b).
In scenario 2 ( Fig. 3e and 7d), up to 115 • C both small-grained magnetite and large high-titanium content grains are demagnetized, both of which carry the pTRM. Between 115 • C and 180 • C, the situation is more complex: large titanium-rich grains carrying the ChRM are demagnetized, together with small-grained lowtitanium magnetite carrying the pTRM as well as intermediate strongly dependent on the grain distribution: they lie on the diag-onal describing the boundary of the titanomagnetite-content/grainvolume distribution. Depending on the distribution the two points may occur at different Curie temperatures and grain-sizes.
Compared to scenario 2 (pTRM), scenario 3 (VRM) shows a similar picture in the small-grained, magnetite-rich half of the diagram. It is observed that less medium-sized and large titaniumrich grains acquired the VRM, compared to the pTRM. This is the reason that the demagnetization plots (Fig. 4) show a significantly weaker VRM than pTRM; the titanomagnetites are less responsive to VRM acquisition than they are to TRM acquisition. The same effect occurs for the broad grain distribution but the effect is considerably more smeared out due to the distribution shape.

AF demagnetization
The effects observed during thermal demagnetization are even more pronounced in the AF demagnetization data: first, as AF demagnetization is done on a timescale of 10 ms (at 50 Hz), the difference in time between acquisition and demagnetization is even larger than in thermal demagnetization, and second, the shape of the AF demagnetization curves is given by a different equation (eq. (11)); the slope of the dotted lines in Fig. 7 indicating the AF blocking condition for different peak fields H is shallower than the dashed lines for thermal demagnetization. Increasing the peak AF field has a similar effect as increasing the temperature: shifting the lines to the top right while rotating them anti-clockwise: titanomagnetites are more responsive to both thermal demagnetization (increases in temperature) and to AF demagnetization (increases in peak AF field) than magnetite, which is more responsive to VRM acquisition and decay (increases in time).
As these effects are more pronounced than for thermal demagnetization, four apparent components arise in the case of the bimodal distribution (Fig. 3): An example is shown for scenario 1, where up to 25 mT, mostly titanomagnetite pTRM is demagnetized (Fig. 7a), and up to 55 mT, mostly titanomagnetite ChRM is demagnetized. In the latter range, small-grained magnetite is also demagnetized, causing an overlap between the two components in the Zijderveld plot. As these grains are small, however, their magnetic moment is weak, and the direction of the pTRM dominates. Above 55 mT, most titanomagnetite grains have no remanence left, and the demagnetization of small-grained magnetite carrying the pTRM dominates. Above 75 mT, the larger-grained magnetite carrying the ChRM are demagnetized.
Similarly for the continuous grain distribution, the two minerals (high-titanium and low-titanium titanomagnetite) show completely separate components in the demagnetization plots (Fig. 4). In addition to plots showing four apparent components ( Fig. 4b  and 4h) similar to the bimodal grain distribution, a further effect is encountered in scenario 2 (Fig. 7d): Here the two components completely overlap from 23 mT, but in different proportions. While further increasing the peak AF, the grain population is progressively demagnetized from two sides: large-grained, high-titanium and small-grained, low-titanium. While the small-grained, lowtitanium magnetite carries only the pTRM, the large-grained, hightitanium titanomagnetite carries both the pTRM and the ChRM. This situation continues up to 72 mT, where both components are simultaneously completely demagnetized. Therefore, the AF demagnetization plot (Fig. 4e) appears to show three components, a low-coercivity component in the pTRM direction, an intermediate coercivity component in the ChRM direction, and a high coercivity component in an intermediate direction. This interpretation is, however, incorrect, as (1) the apparent intermediate coercivity component results from the fact that in this coercivity range the magnetic moment of the large-grained titanium-rich titanomagnetite is relatively larger than that of the small-grained titaniumpoor magnetite, making the direction appear close to the ChRM- Fig. 8. Plot of the cooling rate effect for linear cooling in Thellier and Thellier (1959) type experiments for the three grain distribution cases and model validation case of pure magnetite. Linear fit lines are indicated. direction, and (2) the high coercivity component results from an overlap of medium-sized titanium-rich grains carrying the ChRM and medium-sized iron-rich magnetite grains carrying the pTRM, that both have similar sizes and hence similar magnetic moments, yielding an intermediate direction. Both apparent directional components therefore strongly depend on the grain distribution: the directions obtained from such Zijderveld plots are equally dependent on both the grain distribution and on the directions of the magnetizing fields. Given this particular grain distribution, it becomes obvious that the highest coercivities need not represent the ChRM: in this case the intermediate coercivities are closer to the real ChRM direction due to less overlap, whereas the high coercivities are far less useful for paleodirection reconstructions.

Paleointensity
The Thellier andThellier (1959)/Coe (1967) paleointensity simulations ( Fig. 5 and 6) showed very little dependence on the grain distribution compared to the Zijderveld plots. They are, however, very sensitive to the cooling rate / acquisition time of the remanence, in accordance with previous studies (Halgedahl et al., 1980;Fox and Aitken, 1980;McClelland-Brown, 1984;Bowles et al., 2005;Biggin et al., 2013); lower cooling rates and hence longer acquisition times lead to a steeper slope in the Arai plots. Average slopes have been calculated for a single TRM acquired at various cooling rates from 10 min to 1 Ma (Fig. 8). The relationship between cooling rate (linear cooling) and Arai plot slope is approximately linear (on a logarithmic scale), but varies slightly depending on the grain distribution. For cooling over a day, an NRM/pTRM ratio of ∼1.12 is obtained and for cooling over 1 Ma a ratio of 1.5 to 1.6 depending on the grain distribution. These values correspond to an overestimate of paleointensities of 12% and 50-60%, respectively, without appropriate cooling rate corrections. These values are similar to those obtained by Halgedahl et al. (1980) analytically for slowly cooling rocks: Halgedahl et al. (1980) calculated pTRM acquisition values that would lead to 10% overestimates of paleointensities for a rock cooled over 2 days and about 45% for cooling over 1.6 Ma; the first value coinciding closely with the one obtained here, and the second value being slightly lower.
The slope of the Arai plots is determined by the relative strength of NRM loss during the first heating cycle versus pTRM gain during the second heating cycle. Both the relative independence from mineralogy and the strong dependence on cooling rate can be explained this way: In standard Thellier and Thellier (1959)/Coe (1967) experiments, both heating cycles have the same heating and cooling rates and the same hold-time. Therefore, both cycles activate the same set of grains. The slope of the Arai plot is therefore largely independent of the grain distribution. On the other hand, the total NRM moment that is carried by the set of grains that are activated during the first heating cycle depends on the acquisition time: longer acquisition times generally lead to a higher remanence carried by the magnetic grains and therefore lead to a steeper Arai plot. Minor variations of the Arai plots with mineralogy are due to the fact that the different magnetic minerals have different levels of responsiveness to this cooling rate effect: magnetite is able to continue to acquire a thermoviscous remanence at temperatures below its blocking temperature, thereby increasing its magnetic moment upon slow cooling. Larger grained high-titanium titanomagnetite, on the other hand, block close their Curie temperatures and do not significantly increase in magnetization upon slow cooling. Therefore the cooling rate effect is slightly stronger in titanium-poor minerals, in which case a slightly larger correction factor must be applied (Fig. 8).
The simulations suggest that for slowly cooling rocks a cooling rate correction of up to a factor of 1.5 (for broad grain distributions) to 1.6 (for pure magnetite) may need to be applied. The correction factors obtained here (Fig. 8) agree well with those obtained by Halgedahl et al. (1980). These theoretical predictions have been experimentally confirmed for SD samples by various studies: 6%-12% overestimates for archaeological baked clays refired and cooled over 7 h in the laboratory (Fox and Aitken, 1980); 15% overestimates for synthetic SD magnetite with NRM acquisition on cooling 50 times slower than the Thellier andThellier (1959) experiments (McClelland-Brown, 1984) (equivalent to 8 h in Fig. 8); 11-26% overestimates for remelted volcanic glass containing SD magnetite on 75-fold lower NRM acquisition cooling rate (Ferk et al., 2010) (equivalent to 12 h in Fig. 8); 5-10% overestimates for SD low-Ti titanomagnetite volcanic glasses at 34-fold lower NRM acquisition cooling rate (Leonhardt et al., 2006) (equivalent to 6 h in Fig. 8). Similar values were obtained by others for baked clays and volcanic glasses in the SD range (Papusoi, 1972;Chauvin et al., 2000;Bowles et al., 2005;Yu and Tauxe, 2006) and the PSD range (Yang et al., 1993;Biquand, 1994;Genevey and Gallet, 2002;Genevey et al., 2003;Morales et al., 2006). Biggin et al. (2013) found that the cooling rate effect is weaker for interacting SD, PSD and MD grains than for non-interacting SD grains, with a ∼3% increase in TRM magnitude per order-of-magnitude decrease in cooling rate. In summary, the cooling rate effect on paleointensities in this study coincides well both with theoretical predictions by Halgedahl et al. (1980) and with experimental observations and is more important than the grain size/composition distribution.

Conclusions
The simulations have shown that the presence of mixtures of titanomagnetites has very significant effects on the vector demagnetization plots in all cases except the one where the demagnetization timescale is equal to the acquisition timescale. In particular, two cases can be observed in stepwise thermal demagnetization, one that shows an apparent third component at intermediate temperatures that arises from an overlap of a remagnetization carried by small-grained iron-rich magnetite and large-grained titaniumrich titanomagnetite, and one that shows a continuous curvature between the two components. In both cases, the blocking temperatures of the "intermediate component" are a function of the grain distribution, the acquisition time and temperature and the demagnetization time. In particular, although in clearly bimodal distributions, where two clear distinct Curie temperatures can be measured in the M s (T ) curves, the upper and lower blocking temperatures can be attributed to the two grain populations with two distinct Curie temperatures, in more continuous distributions the exact mineral (Curie temperature) and grain size that the upper and lower blocking temperatures of the apparent intermediate component correspond to is not easily determined. Instead, it depends on the shape of the grain distribution, with the blocking temperatures corresponding to neither of the end-members of the distribution.
The simulated AF demagnetization experiments show particularly strong deviations from the case of a unique magnetic mineral, which is due to the short timescales involved in AF demagnetization and to the different blocking mechanism. For remanent magnetizations with long acquisition times, Zijderveld plots of AF demagnetization experiments may show three to four components. The highest coercivity component is not necessarily equivalent to the primary remanence and does not necessarily correspond to the highest temperature component in an analogous thermal demagnetization experiment. Although the interpretation of such Zijderveld plots is not straightforward, the magnetic remanences carried by different magnetic minerals may appear completely separate in AF demagnetization, which may allow to isolate the paleomagnetic directions of interest.
For paleointensity experiments it was found that the grain distribution affects the slope of Arai plots, but is negligible compared to the effect of the cooling rate of NRM acquisition. The simulations suggest that for slowly cooling rocks a cooling rate correction of up to a factor of 1.5 (for broad grain distributions) to 1.6 (for pure magnetite) may need to be applied. It was also shown that VRM acquisition impacts Arai plots, even though their direction may be indistinguishable from the ChRM. Contrary to directional analysis, paleointensities can be relatively easily analyzed using the cooling rate / Arai plot slope correction factors in Fig. 8. The cooling effect in our simulations is similar in magnitude as theoretically predicted by Halgedahl et al. (1980) and consistent with experimental observations. All this shows that it is critical to identify the presence of mixtures of different magnetic minerals when interpreting demagnetization data for paleomagnetic field reconstruction. Although the same information about the magnetic history of a sample is preserved in mixtures as in pure materials, its interpretation is significantly complicated. Mixtures of different minerals can often be identified from M s (T ) curves: all curves in Fig. 2 significantly deviate from pure magnetite or single-titanomagnetite curves, either showing more than one clearly distinguishable Curie temperature or showing a strong decay at low temperatures leveling off at high temperatures. When such mixtures are identified, additional information about the acquisition times is needed to correctly identify primary magnetizations in Zijderveld plots.