Was the Devonian geomagnetic field dipolar or multipolar? Palaeointensity studies of Devonian igneous rocks from the Minusa Basin (Siberia) and the Kola Peninsula dykes, Russia

V.V. Shcherbakova,1 A.J. Biggin,2 R.V. Veselovskiy,3,4 A.V. Shatsillo,3 L.M.A. Hawkins,2 V.P. Shcherbakov1,3,5 and G.V. Zhidkov1 1Geophysical Observatory Borok IPE RAS, Russia. E-mail: shcherbakovv@list.ru 2Geomagnetism Laboratory, Dept. Earth, Ocean and Ecological Sciences, University of Liverpool, United Kingdom 3Institute of the Physics of the Earth RAS, Russia 4Geological Department, Lomonosov State University, Moscow, Russia 5Institute of Geology and Petroleum Technologies, Kazan (Volga region) Federal University, Russia


I N T RO D U C T I O N
The Earth's magnetic field is known to vary its behaviour considerably on timescales of tens to hundreds of million years, changing its reversal frequency and its average strength (e.g. Biggin et al. 2012;Tauxe & Yamazaki 2015;Meert et al. 2016). The causes and limits of these variations are poorly understood but are key to understanding the dynamics of the deep Earth and the history of geomagnetic shielding of the atmosphere and life from deleterious solar wind effects.
Previous palaeomagnetic studies of Devonian rocks have suggested a complex and unstable Devonian palaeofield configuration.
The first measurements from the Lower Devonian 'Old Red Sandstone' (ORS) in Scotland by Creer & Embleton (1967) showed that two rather different poles (DI and DII) can be isolated for this period. To explain this, it was hypothesized that the DI direction represented a Permo-Carboniferous re-magnetization while the DII direction was considered as true Devonian pole as it agrees with the palaeoclimate zones for the East European platform. Besides DI and DII, other quasi-stable palaeofield directions of normal or reverse polarities interpreted as transitional to DII (Creer 1968;Khramov et al. 1974) were also found in ORS sections.
Later investigations on Devonian rocks from North America and East-European platform suggested that the problem of obtaining C The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. 1265 1266 V.V. Shcherbakova et al. reliable Devonian poles is global (Sallomy & Piper 1973;Khramov et al. 1974;Latham & Briden 1975;Kono 1979;Roy et al. 1983;Torsvik 1985;Jackson et al. 1988;McCabe & Elmore 1989;Orlova 1992;Smethurst & Khramov 1992).
Devonian rocks from Siberia and other parts of Asia have been intensively investigated by a number of researchers over the last 15 yr (Kravchinsky et al. 2002;Bazhenov & Levashova 2011;Orlov & Shatsillo 2011;Veselovskiy & Arzamastsev 2011;Konstantinov & Stegnitskii 2012;Bazhenov et al. 2013;Veselovskiy et al. 2013Veselovskiy et al. , 2016Shatsillo et al. 2014). Most of the objects studied in these papers revealed Devonian palaeodirections but they also suggested the possibility of unusual geomagnetic behaviour in the Devonian.
Given these difficulties in describing Devonian palaeofield directions, it would be extremely useful to compare them with results of palaeointensity (B anc ) studies for this period. However, current palaeointensity databases (wwwbrk.adm.yar.ru/palmag/database e.html; http://earth.liv.ac.uk/ IAGA/) contain only seven published data for the Devonian (Briden 1966a, b;Smith 1967;Schwarz & Symons 1969;Kono 1979;Didenko & Pecherskiy 1989;Solodovnikov 1996). These studies were carried out on different rocks (gabbro, welded tuff and argillite, red porphyritic lava flows, etc.) from sections sampled in Asia, North America and Europe and four different methods were used for these investigations: Shaw, van Zijl, Wilson and Thellier-Thellier. Remarkably, the majority of these data point to the domination of low and very low intensities (see Section 6) which fits with the idea of a peculiar Devonian field inferred from the palaeodirectional data. However, palaeointensity data are too sparse and fall short of current standards of reliability: only data from Solodovnikov (1996) were obtained using the Thellier method but these experiments did not include pTRM checks. Nevertheless, the consistently low results obtained over a variety of regions and lithologies supports that palaeointensity data also indicate that the Devonian field configuration was indeed quite unusual.
It is clear that in order to obtain more definite conclusions on this subject, new Devonian palaeomagnetic and palaeointensity data that satisfy modern criteria of reliability are badly needed. Here new palaeointensity results from rocks of this age sampled in the Minusa Basin, Siberia and the Kola Peninsula are presented and discussed. We will first outline our general sampling and experimental results before presenting geological backgrounds and results from each of the studied regions in turn.

S A M P L E C O L L E C T I O N A N D E X P E R I M E N TA L P RO C E D U R E S
Rocks from the Minusa Basin (Siberia) and the Kola Peninsula were sampled by groups led by Dr Shatsillo andDr Veselovskiy, respectively, from 2007 to 2015. Although the palaeomagnetic directional studies are already reported elsewhere (Shatsillo 2013;Veselovskiy et al. 2013Veselovskiy et al. , 2016, we briefly describe the experimental procedure by which they were obtained since their results play an important role in the further discussion. The procedures of sampling and measurement were very similar. Samples were taken using a portable drill (Pomeroy Ltd.) and oriented using sun and magnetic compasses. In the Kola Peninsula (Veselovskiy et al. 2013(Veselovskiy et al. , 2016, up to 24 oriented samples were taken from dykes with thicknesses of between 1 and 60 m and 6-12 samples from dykes with thickness <1 m. In the Minusa Basin, 6-10 samples were taken from each sampling site. The requirement to take up to 24 samples from a thick dyke aims to get a representative palaeomagnetic record from different parts of the dyke (Van der Voo 1990). However, as a dyke of 0.5 m thickness is much more homogeneous in composition the number of collected samples can be reduced. To avoid the effects of lightning strikes and unaccounted rotations of lava blocks, in each lava flow samples were collected over an area of several square meters and/or were laterally spread in a given lava flow over a distance of up to 5-20 m.
Palaeomagnetic directional studies of Kola samples were made in Moscow at the Lomonosov Moscow State University, the Geological Institute of the Russian Academy of Sciences (RAS), and the Institute of the Physics of the Earth (IPE RAS), and also at MIT (USA). The palaeomagnetism of rocks from the Minusa Basin was investigated in the IPE RAS. The main method to isolate the characteristic remanent magnetization (ChRM) was stepwise thermal demagnetization performed with not less than 12 temperature steps up to 580-700 • C using the nonmagnetic furnaces MMTD-80 and ASC TD-48. For comparison, a few experiments of alternative field's demagnetization (AF-demagnetization) were undertaken on sister samples. For these directional analyses, the remanent magnetization was measured with an AGICO JR-6 spinner magnetometer and/or with a 2G Enterprises SQUID magnetometer.
All rock magnetic and thermal palaeointensity investigations were carried out in the Borok Observatory of IPE RAS with nonoriented sister cubes of 1 cm in edge length. The Curie points T c and the thermal stability of magnetic minerals were estimated from thermomagnetic curves of strong-field magnetization M s during heating-cooling cycles to incrementally higher temperatures T i . These measurements were made with a Curie balance in an external magnetic field B = 450 mT. Hysteresis loops and additional high-field thermomagnetic curves (up to 900 mT) were measured using a variable field translation balance (VFTB) produced by Petersen Instruments company. To estimate the magnetic hardness and mineralogy of samples, measurements of magnetic susceptibility, hysteresis loop parameters, such as coercive force B c , remanent coercive force B cr , M s and remanent strong-field magnetization, M rs , were made. The domain structure (DS) was estimated from the Day plot (Day et al. 1977) using ratios M rs /M s and B cr/ B c and also from the thermomagnetic criterion by measuring the tail of partial thermoremanent magnetization (pTRM) surviving the thermal demagnetization of the given pTRM as proposed by Shcherbakova et al. (2000).
Palaeointensities were determined using three different techniques. Palaeointensity determinations by the classical Thellier-Coe method were undertaken at the Borok Observatory of IPE RAS (Thellier & Thellier 1959;Coe 1967). These experiment consisted of a sequence of paired heatings in 1 atm air to a set of increasing temperature T i , i = 1, . . . n. The first heating-cooling step to T i takes place in zero field, the second heating is also performed in nonmagnetic space followed by cooling in the laboratory field, B lab . This value was set to 20 µT for the majority of the experiments but values of 30, 5, 3 or 1 µT were also used for verification purposes. Double heatings were carried out over at least 15 steps up to 650-690 • C; pTRM checks and room-temperature susceptibility measurements were performed after every second step. The results of the Thellier procedure were used for constructing the Arai-Nagata and Zijderveld plots (in specimen coordinates), which enabled us to control the choice of the appropriate temperature interval (T 1 , T 2 ) for estimating B anc in the subsequent analysis. The temperature range (Т 1 , Т 2 ) for estimating B anc from the Arai-Nagata plots was selected as close as possible to the interval where the ChRM was isolated in oriented sister-samples analysed in the palaeodirectional experiments.
Two to four sister specimens were used from each sample for palaeointensity analyses. Some of them were subjected to the Thellier-Coe procedure in a full-vector three component vibrating sample magnetometer (3D-VSM) or in a two component vibrating sample magnetometer (2D-VSM), both constructed at the Geophysical Observatory 'Borok'.
The sensitivity of the 3D-VSM is 10 -8 Am 2 and the maximum available external field is 0.2 mT. The 2D-VSM allows demagnetization and remagnetization of (p)TRMs while measuring continuously the horizontal components of specimen's remanent magnetic moment. The noise threshold of the 2D-VSM is 3 × 10 −9 Am 2 and the maximum available external field is 0.2 mT (residual field < 100nT). Specimens not treated in a VSM were heated in an electric furnace with a residual field <50 nT and were measured with the JR6 magnetometer. No systematic variations of the palaeointensity result were observed between devices.
The 2-D VSM was used for haematite samples where the magnetization was too weak for the 3-D device. In parallel, sister-cubes of these samples were also measured on the JR6 magnetometer. To minimize the consequences of measuring only two components, specimens were placed in the holder so that these constituted the two strongest components. We selected only those specimens whose weakest (vertical and unmeasured) component Z was at least 10 times less than the total intensity of the 2-D horizontal vector √ (X 2 + Y 2 ). From here one can see that the error in intensity is Z/2 √ (X 2 + Y 2 ) which amounts to not more than 5 per cent. Palaeointensity results for sister-cubes from the same sample from the 2-D and JR6 magnetometers were indistinguishable.
In addition to the Thellier-type experiments, Wilson's method (Wilson 1961) of palaeointensity determinations was also applied in the Borok Observatory. In this procedure, thermal demagnetization curves NRM(T, B anc ) and full thermal remagnetization curves TRM(T, B lab ) are compared to find the temperature interval, (Т w1 , Т w2 ), (Т w1 < Т w2 ), where the shape of the two curves is the most similar. To perform the comparison, we introduced a fitting function TRM * (T) with a help of a multiplying coefficient k w defining TRM * (T) = k w × TRM(T). If such an interval (Т w1 , Т w2 ) is found, the palaeointensity estimate (B anc * ) can be obtained from the coefficient of similarity, k w = NRM(T)/TRM(Т) = B anc * /B lab , T ∈ (T w1 , T w2 ). We highlight that the similarity of NRM(T) and TRM(T) curves provides a strong argument in favour of the thermoremanent nature of the NRM. The criteria proposed by Muxworthy (2010) were used to judge whether a Wilson palaeointensity experiment was successful or not.
The third method used was the microwave IZZI -Thellier protocol (Tauxe & Staudigel 2004) performed using the third generation microwave system (Walton et al. 1996;Hill & Shaw 1999 at the University of Liverpool. This system uses high frequency (14 GHz) microwaves to directly stimulate the magnetic grains in a 5 mm cylindrical rock sample, simulating thermal demagnetization/remagnetization while minimizing the thermal energy within the bulk sample. The magnetic moment of the sample is measured using a Tristan 3-axis SQUID magnetometer mounted in-line with the resonant cavity. A 3-axis coil assembly around the cavity allows the magnetic field vector to be fully defined in the range 0-0.1 mT. Microwave experiments may be performed with a similar range of  versatility as thermal counterparts, using increases in the microwave power and duration of application in place of temperature, and have been demonstrated many times to give equivalent results (e.g. Biggin et al. 2007;Hill et al. 2008;Poletti et al. 2013;Monster et al. 2015).
The IZZI (in field -zero field -zero field -in field; Tauxe & Staudigel 2004) protocol of the Thellier experiment follows the same principles as Coe's modification (see above) but reverses the order of the demagnetization/remagnetization at each temperature or microwave power-time step. This serves to introduce a characteristic zigzag shape to the Arai-Nagata plot where MD-type effects are present negating the need for pTRM tail checks (Riisager & Riisager 2001). Repeated pTRM-checks for alteration were made at each alternate power-time step.

PA L A E O I N T E N S I T Y S E L E C T I O N C R I T E R I A
Parameters associated with the Thellier-type experiments (both thermal and microwave) were calculated according to the Standard Palaeointensity Definitions ). Our selection criteria take their lead from the SELCRIT2 criteria (Biggin et al. 2007) which were shown to produce an accuracy rate (±10 per cent) of 78 per cent and produce mean deviations of less than 5 per cent accuracy when applied to large simulated datasets from materials subject to MD-like effects ). The microwave system occasionally produces problems with regards to reproducing absorbed power intergrals that can contribute noise to Arai plots and pTRM checks in particular. We therefore decided to relax slightly the DRAT criterion of SELCRIT2 but to compensate for this by introducing a new CDRAT criterion.
1. n, the number of points used to fit the straight line on the Arai-Nagata plot ≥4.
2. FRAC, the NRM fraction used for the best-fit on the Arai-Nagata plot determined entirely by vector difference sum  Enkin (2003). X indicates that there were insufficient data, symbol '?' indicates that the result was inconclusive; percentage of unfolding at maximum accuracy of distribution of the precision parameter k. * -reversal test is positive (see text). 3. β, the standard error on the slope of the best-fitting line normalized by its slope ≤0.1.
6. MAD, the maximum angular deviation of the selected component on the accompanying Zijderveld plot ≤15 • . 7. α, the angle between the origin-anchored and floating fit to the selected component on the accompanying Zijderveld plot (Selkin et al. 2000) ≤15 • (although note that it was not possible to calculate this in cases where only two components of remanence were measured).
8. DRAT, the maximum difference ratio measured from pTRM checks ≤±15 per cent (note that SELCRIT2 used ±10 per cent but our leniency is counterbalanced by incorporation of the next criterion). 9. CDRAT, the cumulative DRAT calculated from summing from all pTRM checks up to the maximum temperature step used to obtain the best-fitting straight line ≤±16 per cent (note that CDRAT was not included in the original SELCRIT2 criteria).
10. DRAT TAIL , the equivalent of DRAT for pTRM tail checks (note that these were only performed in a minority of our experiments) ≤±10 per cent.
The moderate failure rate (45 per cent) experienced in this study resulted in several sites being represented by only one or two specimen result. In such cases, we retained the result for the purpose of discussion but did not calculate a virtual dipole moment (VDM) for incorporation into the global database. Where three or more specimens gave acceptable results, we calculated a site mean palaeointensity and used the inclination associated with the site mean direction to calculate a VDM which we then evaluated using the Q PI criteria of .

Geological background and sampling
Sampling of rocks in the Minusa Basin was performed within the Minusinskaya Valley region. The sections are formed by thick (1-2 km) lava piles of Early and Middle Devonian age which are part of the Byskarskaya series. More than one hundred distinct volcanic and subvolcanic units from four mafic sections Tuba, Koksa, Sisim and Truba of this series were sampled for palaeomagnetic investigations (Fig. 1).
The Minusa Basin is usually treated as a Middle-Late Palaeozoic rift structure imposed on the epi-Caledonian Siberian continent and being contiguous with the Siberian craton since the Middle Palaeozoic. A part of the Byskarskaya series of Minusa Basin is subdivided into a set of formations (Immirskaya Fm., Tomskaya Fm., Chilanskaya Fm., etc.). Lavas from this series have a bimodal composition of rhyolites and alkali-rich basalts (accompanied by tuffs of the same composition). Some of the studied outcrops are located on the wings of the large-scale folds that allow fold tests to constrain the age of magnetization relative to end-Palaeozoic aged folding deformations (Buslov 2011).
There is a set of modern isotopic dates obtained from U-Pb and Ar-Ar methods for separated sections of Byskarskaya series, yielding ages from 407.5 ± 0.2 to 386 ± 4 Ma (Vorontsov et al. 2013).The lavas are covered by sediments whose faunal remains are estimated as Eifelian -Givetian (398-385 Ma in age). Hence, the age of investigated sections should correspond to Emsian-Eifelian (408-388 Ma). The Immirskaya Fm from the Sisim section contains the oldest stratigraphy units studied here.

Palaeomagnetism
Preliminary palaeomagnetic directional results are published by Shatsillo (2013) and will be summarized here. The majority of samples used for palaeointensity estimates produce clean high temperature (300-680 • C) components ( Fig. 2) that yield mutually consistent directions for each site studied (e.g. sites selected for palaeointensity study mostly had N > 5 and k > 50).
Principal component analysis of the results of thermal magnetic treatment reveals a wide spectrum of the NRM components ( Table 1). Two of them, named here as 'N' and 'S' (Fig. 3) and combined into 'green group', have been found in almost all of the studied localities, separated by a distance of more than 200 km. Site-mean directions of these components passed the reversal test (γ /γ cr = 8.8/10.2; McFadden & McElhinny 1990). The fold test (Enkin 2003) is positive for N component and inconclusive for S component due to the similarities of the lavas' bedding. Moreover, the corresponding VGPs match well with the Siberian platform poles reported for the Lower Devonian by Bachtadse et al. (2000). Based on the high quality of the directions at sample and site levels and taking into account the evidence of their primary nature, presented above, we consider that all samples with N and S components of NRM selected for palaeointensity analysis likely record the primary Devonian direction.
Besides the N and S components, up to five ChRM components of NRM can be found in most of the sites ('E', 'NW', 'NWSE', 'SW' and 'P', Table 1). These components occupy a high blocking temperature interval and their NRM vector behaviour during the thermal treatment usually looks similar to those of N and S components. NW, NWSE and P components successfully passed the fold test, so the time of their origin can be estimated as pre-Late Palaeozoic (pre-Hercynian). Corresponding VGPs of all these  -13(f, site11-13) for Kola collection. Full lines are the continuous temperature curves of the pTRMs(T 1 ,T 2 ) after switching off the lab field 100 µT at T 2 and cooling toT 0 , followed heating to T 1 and cooling from T 1 to T 0 ; all are made in zero field. The arrows point the direction of temperature change. Vertical lines mark the temperature interval in which pTRM(T 1 , T 2 ) was created. components (except P-component) are offset from the known Phanerozoic palaeomagnetic poles of Siberia . Thus the discussed ChRM components could not be related with later remagnetization. Similarly, they are unlikely to be due to unaccounted tectonic rotations because they are consistent over a large area and exist at the same outcrops (but not at the same lava flows) where the primary N and S components are. There are therefore convincing arguments that the E, NW, NWSE and SW components of NRM are of primary Devonian age. We therefore use the corresponding samples for palaeointensity determination, but with caution and combine them in a 'yellow group'.
The last 'red group' of ChRM components can be distinguished in studied samples from Sisim and South Koksa sections and include directions that appear as outliers (Table 1)   direction while showing high levels of internal consistency (N = 6, k = 64-287). These components have been found in high temperature interval (450-700 • C). The direction of these ChRM components is stable for all samples within each site, but site-mean directions have a quasi-chaotic distribution on a stereoplot and each direction from this group can be met in one site only. They have been found at the same outcrops (but not in the same flows) where we observed the primary components N, S, NW and SW. They are unlikely to be a result of the later remagnetization event and, possibly, similar to the NW, SW, NWSE and E components have a primary origin that reflects an anomalous direction of Devonian palaeofield. We used the samples with 'outlier' components for palaeointensity estimation, but treat their results with great caution.  The weak Devonian magnetic field 1275   Besides these sites, one more NRM component, named 'P', is also included in the 'red group' because its pole lies close to a known Jurassic palaeomagnetic pole of Siberia ) and therefore may be contaminated by a later remagnetization event.

Rock magnetic properties
Experiments measuring the temperature dependence of strong-field magnetization (M s ) were performed on a sister-specimen of each sample measured for palaeointensity. For the majority of samples, M s is weak ranging from 10 to 200 A m -1 . Sets of successive thermomagnetization curves M si (T) of samples from the Minusa collection fall into three basic types (Figs 4a1-c1). Samples from the Sisim sites 041 and 044, and Tuba site 030 show thermomagnetic curves typical for haematite (Fig. 4a1) and have the lowest intensity. The haematite presence is confirmed also by their very high magnetic hardness so that even fields up to 1 T cannot saturate these samples (Fig. 4a2) which makes it impossible to define the hysteresis parameters M rs /M s and B cr /B c. Samples from Sisim sites 044 and 045, Truba sites 031 and 032, Tuba site 028, and Koksa-S sites 007 and 009 also demonstrate thermally stable curves M si( T), but in addition to haematite, the presence of magnetite is also prominent (Fig. 4b1). These conclusions agree with the wasp wasted (though not yet fully closed) shape of hysteresis loops for these samples suggesting a mixture of magnetically soft (magnetite) and hard (haematite) contributions to the magnetization (Figs 4b2 and c2). Samples from Sisim site 053, Truba site 33, and Koksa-N sites 70-1, 70-5, 70-6 and 70-7 also showed a mixed magnetomineralogy and some degree of thermal instability to heating (Fig. 4c1). Successful palaeointensity results were obtained from sites associated with all three of the thermomagnetic curve types.
The presence of magnetite and haematite in Sisim sites 41 and 45 is confirmed by the electron microscopic images of thin sections taken from samples 417 (site 041; Fig. 5a) and 440 (site 45; Figs 5b and c). These figures display a typical structure of high-temperature (heterophase) oxidation of the initial titanomagnetite (TM) characterized by the exsolution of TM into ilmenite lamellae and magnetite cells. The EDS analyses did indeed show that the atomic content of lamellae corresponds to that of ilmenite. However, the results of EDS analysis for cells was not conclusive as the chemical formulas of magnetite and haematite are insufficiently well-spaced. To remove this uncertainty, we carried out etching of thin sections which helps to distinguish haematite and magnetite cells since etching with hydrochloric acid dissolves magnetite but not haematite. As is seen (Fig. 5a) the cells in the grain from sample 417 (site 041) were not dissolved proving that they consist of haematite only as it was suggested by the thermomagnetic analysis (Fig. 4a1).
This fact indicates that the process of oxidation had deeply proceeded in this sample resulting in transformation of magnetite to haematite.
Conversely, sample 440 was not fully oxidized to haematite as some cells were dissolved by etching as it is demonstrated in thermomagnetic criterion (Shcherbakova et al. 2000;Fig. 6). Note that only samples which display the hysteresis loop of the type shown in Fig. 4(c2) are presented on the Day plot as the loops of other samples are far from being saturated. As is seen, the distribution of the representative points on the Day diagram is confined to SD-small PSD regions for magnetite suggesting the prevalence of fine magnetic grains. The same conclusion can be made from plots 6b to 6f as the pTRMs created at temperature intervals where the palaeointensity was estimated have no tails supporting SD-like behaviour. This conclusion agrees well with the electronic microscope observations (Fig. 5) which showed a great abundance of exsolution structures subdividing coarse TM grains into ilmenite lamellae and magnetite cells.

Palaeointensity determinations
Palaeointensity experiments were performed on 109 specimens from the Minusa Basin collection in the Geophysical Observatory 'Borok' and the University of Liverpool. In all, 61 specimens from 15 different sites of the five investigated sections passed the criteria listed above and we will consider these samples only. Fig. 7 provides examples of Arai-Nagata, orthogonal and Wilson's plots for samples representing the cases discussed above of pure haematite and mixtures of magnetite and haematite. For all samples shown in Fig. 7, no significant sample alteration was observed over the entire temperature intervals as indicated by positive pTRM checks. Most samples also have a minor low-temperature secondary component as shown in the orthogonal plots in Fig. 7 and so these temperature intervals were not taken into account.
The six samples from site 041 and four samples from site 044 from the Sisim section differ from other sites since they contain haematite as the main magnetic mineral as illustrated in Figs 4(a1)-(a2) and Fig. 5. Correspondingly, for these samples the palaeointensity was determined for high temperature intervals up to 690 • C and these yielded very low mean B anc values of 8.8 and 5.6 µT (Table S1, and Figs 7a1, c1).
Remanence carriers in samples from other sites are mostly a mixture of magnetite and haematite. The Arai-Nagata diagrams for most of these samples allowed large temperature intervals from 300-400 to 600-685 • C (depending on sample) to be used to calculate the palaeointensity from a combination of both these minerals. Examples of such fits are shown in Figs 7(b1) and (c1). Example of successful Wilson experiment results are also given in Figs 7(a3)-(d3) while Figs 7(e1) and (e2) presents examples of successful microwave Thellier results.
Nearly half of the samples measured did not produce an acceptable palaeointensity measurement. An example of such a sample is presented in Fig. 8. The formal reason for its rejection is that it does not pass the selection criteria (Table S1). Another obvious reason is that the Arai-Nagata plot of NRM for this sample may be divided in three quasilinear sections (Fig. 8a) having at the same time practically one component NRM vector (Fig. 8b). Finally, this sample exhibits a significant drop in M s (T) intensity (Fig. 8c) after heating to 400 • C that most likely reflects the presence of maghemite or cation-deficient magnetite which on laboratory heating experienced transformation to haematite. This drop in M s (T) likely reflected the drop on the low-temperature part of the Arai-Nagata plot (Fig. 8a).
A summary of results on palaeointensity determinations obtained by both the thermal Thellier (TT) and microwave (MW) methods for the Minusa Basin is given in the upper section of Table S1 and includes 39 successful estimates from all four sections. A further 19 successful estimates produced by the Wilson method are given in Table 2. The order of the sites within given sections reflects their relative stratigraphic positions but, because of problems of regional correlation, the relative ages of the sections themselves, together spanning 15 Ma, are not known.
Eleven site mean intensities of the palaeofield were produced and eight of these were averaged over 3 or more single determinations. All but one of these determinations indicate either low, from 10 to 20 µT, or very low <10 µT field intensity in the Devonian. Remarkably, low intensities were obtained from all four sections (Koksa, Sisim, Truba and Tuba) and from different directional groups (S, N, SW and NWSE) as well as from sites producing outlier directions. The only exception is site 30 (NWSE, 'yellow' component) of Tuba section with the mean value of B anc = 38.2 µT which is closer to the recent field intensity in this region. The Arai-Nagata diagram, orthogonal and Wilson's plots for the sample 125 from this site are shown in Figs 7(d1)-(d3). Another flow (28, 'yellow' component) from this section has the same NWSE direction but low mean intensity B anc = 6.9 µT in accordance with the results obtained for most of sites presented in the upper part of Table S1. Both of these mean palaeointensities included results from all three palaeointensity methods.

Geology, sampling and rock magnetic properties
The Kola Devonian Alkaline Province is located on the NE part of Fennoscandia and includes many dyke swarms and alkalineultramafic plutons such as the Khibiny and Lovozero Massifs. All of these intrusions are hosted in highly metamorphosed rocks of Archean-Proterozoic age. According to a number of Ar/Ar and Rb-Sr isotopic ages (Arzamastsev et al. 2009(Arzamastsev et al. , 2017Veselovskiy et al. 2013), the Devonian magmatic activity over the Kola Peninsula took place between 390 and 360 Ma. Detailed palaeomagnetic studies (Veselovskiy & Arzamastsev 2011;Veselovskiy et al. 2013Veselovskiy et al. , 2016 were performed on more than 2000 samples representing 130 dolerite and alkaline dykes of Devonian age. Some of these showed thermally stable components of magnetization (D) which is bipolar and exhibits low positive and negative inclinations and ENE or WSW declinations. This D component was found in 25 per cent of studied dykes, and its palaeomagnetic pole (Veselovskiy et al. 2013(Veselovskiy et al. , 2016 lies closely to the Early Devonian segment of the Baltica's apparent polar wander path (APWP; Torsvik et al. 2012). Hence, it is likely that the D component represents the primary magnetization of the Devonian dykes. Its primary nature is additionally supported by: (1) positive baked contact test (Fig. 10e)   More than 70 non-oriented samples were passed to the Geophysical Observatory 'Borok' and to the University of Liverpool for palaeointensity experiments. Here we discuss only six dykes (Fig. 9), all characterized by the D component, for which some palaeointensity results were obtained (parts 2 of Table S1 and S2).
Sites 18-13 and 11-15 are from two further alkaline dolerite dykes ( Fig. 9) with thicknesses of 80 and 5 m, respectively. The Devonian age of these dykes is assumed on the basis of the similarity of their strikes and petrographic characteristics to neighbouring dykes, located 1 km to the NW and SE, whose Ar/Ar isotopic dating, obtained from feldspar, is 382.8 ± 5.6 and 378.0 ± 5.6 Ma, respectively (Arzamastsev et al. 2017). The three other dykes (sites 22-09, 28-09 and 25-11, Fig. 9) are 1.5-m-thick alkaline lamprophyres and are located in the southern part of the Kola Peninsula belonging to the Kandalaksha dyke swarm. These have similar petrographic characteristics, and, 50 meters from site 22-09 there is a dyke with a date 375.1 ± of 3.9 Myr determined by the Ar/Ar method (amphibole; Veselovskiy et al. 2013).
Examination of polished sections of samples from this collection revealed large TM grains exhibiting typical magnetiteilmenite exsolution structure resulting from high-temperature oxidation (Fig. 11). To refine the content of these two phases, X-ray powder diffraction (XRD) was performed. This diffraction pattern clearly indicates a presence of two phases: pure magnetite with a cell parameter of a = 8.398Å and a rhombohedral ilmenite phase with a = 5.087Å. The lines are well shaped with no evidence of another spinel phase or significant line broadening. From here, it follows that the matrix consists of fine magnetite cells intermixed with hemoilmenite lamellae. As is seen from this figure, typical sizes of magnetite cells vary from submicron in size to 5 µm and fall in pseudo-single-domain (PSD) range. This agrees with the hysteresis measurements (Day plot, Fig. 6a) that produced M rs /M s values in the range 0.1-0.2 and results of the thermoremanent criteria (Fig. 6f). Thermomagnetic curves of M si displayed classical magnetite shapes and insignificant thermal instability as illustrated in Fig. 12(a2)-(b2).

Palaeointensity determinations
A total of 56 palaeointensity experiments were performed on samples from the Kola collection; 29 of these performed using the TT and MW methods and 17 of these performed using the Wilson method) were successful (parts 2 of Tables S1 and S2). Examples of thermal palaeointensity experiments are shown in Fig. 12 for samples 099 (site 11-15) and 225 (site 18-13) together with orthogonal (a4-b4) and  plots; examples of successful microwave palaeointensity experiments are provided in Figs 12(c1) and (c2). The vectors of NRM within most samples were composed of two or even three low and high temperature components as is seen from the corresponding orthogonal plots. However, the hightemperature component still covers a substantial part of the entire Arai-Nagata diagram allowing a palaeointensity estimation to be made.

Origin of remanence
The data presented here from both Minusa and Kola collections assume that the NRM of the studied samples is of thermal origin but it has long been discussed that the characteristic remanence of igneous rocks may have a chemical or thermochemical origin (CRM or TCRM instead of TRM; Smirnov & Tarduno 2005). This presupposes that during primary cooling or re-heating, mineralogical alterations occurred at temperatures below the Curie temperature of the remanence carrying minerals. The difficulty of identifying the origin of the NRM is aggravated by the fact that they may produce no significant deviations from the linearity of the Arai-Nagata plot that is observed (Smirnov & Tarduno 2005). Draeger et al. (2006) carried out extensive experiments with artificial CRM and TCRM and indeed found that most of Arai-Nagata plots were 'linear over almost the entire TRM blocking temperature range, whether pTRM checks are positive or not'. These conclusions were later confirmed by Gribov et al. (2017). Importantly, these experiments demonstrated that an apparent palaeointensity inferred from the (CRM, pTRM) diagrams is usually underestimated by a factor of approximately two, though the exact bias to the determination strongly depends on how the CRM was imparted and ranges from 20 to 200 per cent.
From the above it follows that only detailed magnetomineralogical studies can provide evidence in favour of thermoremanent nature of NRM. Strong field thermomagnetic curves M si (T) (Figs 7 and 12) suggest vigorous, high temperature oxidation took place in the rocks. Direct electron-microscopic observations (Figs 5 and 11) showed that this oxidation resulted in the high-temperature exsolution structures with micron sized lamellae and cells. Gapeev & Tselmovich (1983, 1986 demonstrated experimentally that lower temperature oxidation tends to produce finer lamellae because the diffusion rate strongly decreases with falling temperature. In particular, they showed that exsolution lamellae with sizes 1 µm and more provides evidence for oxyexsolution of primary TM formed during the initial cooling of subaerial basalts between 700 and 900 • C. This suggests a thermal origin of the NRM and supports the reliability of the results reported here. Although the microscopic analysis was restricted to three sites, rock magnetic analyses were much more widespread (they were performed for each sample) and strongly suggested that these were representative of all studied samples.
Another argument for the thermal origin of remanence is that the same low palaeointensity results were obtained on different sites with different mineralogies. Specifically, as evident from analysis of the temperature intervals of estimation of B anc presented in Table S1 and in Figs 7 and 12, rocks from the most reliable Minusa 'green' group consisting of Koksa-N and Truba sections plus the site 45 from the Sisim section, contain a mixture of (possibly cationdeficient) magnetite and haematite NRM carriers. At the same time, as it was demonstrated in the Section 5.1, carriers of NRM of samples from the Kola collection are presented by magnetite exclusively. Nevertheless, both collections show very similar low values of B anc .

Comparison of methods and reliability of results
A summary of site means B anc for Minusa and Kola collections for each method used here (TT, MW, Wilson), is given in Table 3 and Fig. 13. For various reasons (in particular, because of limited sample material) not all sites produced acceptable results from all three methods. Nevertheless, the following conclusions can be deduced.
1. Though an overall coherence between methods can be observed, there are small (average differences ≤4 µT) systematic offsets between all three of them in terms of the palaeointensity results they produce.
2. The microwave method produced results that were systematically higher than both the thermal Thellier and Wilson methods from samples from the same site. The Wilson method, in turn produced marginally higher estimates than those of the thermal Thellier. These properties are clearly illustrated in Figs 13(a)-(d).
3. The overall success rate for microwave palaeointensity experiments (31 per cent) was somewhat lower than that for the thermal Thellier (50 per cent) and Wilson experiments (81 per cent). This is, in part, due to the fact that microwaves cannot effectively demagnetize or remagnetize haematite.
The high apparent success rate of Wilson experiments is partly due to the pre-selection of the samples made simply by eye: in the cases when the primary and secondary thermomagnetic curves were obviously different, the samples were not analysed and are therefore not even considered as a rejected sample. Both Figs 7(a3)-(d3) and 12(a5)-(b5) exemplify the traditional way of presentation results of Wilson's experiments which provides a direct impression of the similarity of the thermomagnetic curves. Besides, the very shape of the primary and secondary thermomagnetic curves is also quite informative to estimate the presence of different magnetic phases and degree of chemical and structural transformations possibly developed during the heating. But it does not give a formal criteria for objective acceptation or rejection of the results. This important part of the study was performed on the basis of estimation of the relative standard error (RSE) of determination, that is the standard error normalized by the field estimate. Muxworthy (2010) who proposed this criteria, took the critical value RSE = 0.02. In order to increase the reliability of the data, we decreased this value down to 0015. Two examples of diagrams illustrating the way of calculation of these criteria are given in Fig. 14. Next, following Muxworthy (2010), the RSE criterion was completed with the well-known Kolmogorov-Smirnov two sample test. All results for the successful Wilson's experiments are listed in Table 2.
The cause of the small but systematic offsets between the microwave and thermal Thellier palaeointensity experiments is far from clear. Biggin (2010) reported a meta-analysis of paired microwave and thermal palaeointensity studies and concluded that the discrepancy was, in general, in the opposite sense to that observed here (i.e. that microwave results tended to be lower than thermal results) and likely caused by multidomain effects being exaggerated in the thermal results because of the specific protocols chosen. It makes good sense that the same offset should not be present here because both techniques apply double-heating protocols and incorporate checks against such non-ideal behaviour. The presence of this new discrepancy, though small, presents a new challenge that is beyond the scope of the current study to explain. The most pertinent finding is that all three methods consistently indicate that the palaeointensity was weak in the vast majority of these Devonian sites. Table 4 summarizes the key data from the investigated collections including the calculated virtual dipole moments (VDMs) and the associated Q PI criteria . The AGE criterion is met by all estimates on account of the isotopic age constraints and internally consistent palaeomagnetic directional behaviour (though the latter are unusual at the formation level for the Minusa sam- Figure 13. (a-c) One-one-comparisons of palaeointensity estimates produced by the three different methods at the site mean level. Note the logarithmic scales.
Dashed lines indicate factor of 2 discrepancies. (d) Palaeointensity estimates produced by the three methods versus site number. Absence of data means no successful paleointensity determinations were obtained for the given site. ples). Only two of the twelve obtained VDMs have N ≥ 5 and s.d./B anc ≤ 25 per cent as required to meet the STAT criterion. The TRM criterion requires reasonable independent evidence to be presented that the component of remanence is thermal in origin. Based on the microscopic analysis results, the detailed discussion above, and the demonstrated similarity between thermal decay and remagnetization curves, all VDMs that were calculated using one or more successful Wilson experiment were judged to have passed the TRM criterion. The ALT criterion is automatically passed by all site mean VDMs on account of pTRM checks being positive under reasonable criteria. Similarly the MD criterion is met by all on account of at least one estimate comprising the mean being produced by the IZZI approach (without showing appreciable zigzagging), positive pTRM tail checks, and/or having a high (>0.6) associated FRAC value. The ACN criterion conflates the three relatively minor considerations of anisotropy of remanence, cooling rate and nonlinear field dependence of TRM. Cooling rate is not considered an issue here on the basis of the remanence carriers being PSD and/or interacting SD and therefore being largely immune to cooling rate effects . Similarly, both the laboratory and inferred palaeofield intensities were generally <20 µT implying that nonlinear effects are also likely to be very small (Selkin et al. 2007). Anisotropy of TRM was gauged by the γ parameter which measures the angle between the applied field and the last pTRM step used to obtain the palaeointensity estimate . Where the majority of estimates had γ < 3 • , the entire ACN criterion is judged to have been passed. All but two of the VDMs were calculated from estimates produced by at least two independent methods and therefore passed the TECH criterion; given the small but systematic (and as yet unexplained) differences presented in Fig. 13, we consider our multi-method approach to be useful in enhancing confidence in the reliability of the mean estimates. All of our sites were single lithologies with no strong differences in unblocking behaviour preventing the LITH criterion from being met. Finally, all raw palaeointensity data have been made available online at pc http://pcwww.liv.ac.uk/∼biggin/Minusa Kola/ allowing the MAG criterion  to be met for all sites. Overall, the Minusa Basin sites produced VDMs with Q PI values ranging from 5 to 8 out of a possible 9 while those from the Kola Peninsula had values of 6 or 7. These include some of the highest values published to date (although it should be noted that the overwhelming majority of Phanerozoic-aged data still await assessment) and demonstrate that these criteria provide a useful mechanism both to quantify relative reliability, and to motivate the meeting of widely acknowledged standards in palaeointensity determination.

Implications of results
The data presented here, together with all previously reported data on palaeointensity and palaeodirections in the Devonian, suggests a complex geomagnetic field structure, weak in strength and having different unstable VGP positions. Published VDMs for the Carboniferous, Devonian and Silurian from the World Database on palaeointensity are summarized together with our new determinations in Fig. 15; all previous data determined by different methods are shown.
Our new evidence for a weak time-averaged field in the Devonian, together with the recognized tendency of palaeomagnetic directions to be somewhat erratic in this period, is suggestive of a long-term departure from the uniformitarian view of a dipole-dominated field geometry through nearly all of Earth's history (Evans & Hoye 2005). The field is considered as dipole-dominated if the coefficient of bipolarity f dip which determines the relative strength of the dipole field at the outer core-mantle boundary exceeds 0.4 (Driscoll & Olson 2009;Christensen 2010) otherwise it is defined as having a multipolar geometry. As shown by Aubert & Wicht (2004), Morin & Dormy (2009), Driscoll (2016 this scenario can be realized due to a decrease in convective power which promotes a transition from a dipolar strong-field state to a non-axial, weak-field dynamo state with significant multipolar components. Conversely, geodynamo simulations have also demonstrated that it is possible to achieve a similar transition in the behaviour by increasing the convective power such that the dynamo enters a strong-field multipolar state (Olson 2006;Driscoll 2016). A similar effect can also apparently be achieved without changing the net forcing at all but rather increasing the heterogeneity of core-mantle forcing or preferentially extracting heat from the core in the equatorial regions (Olson et al. 2010).
Dynamo predictions of behaviour produced by both weak-field and strong-field multipolar simulations broadly agree with the struc-ture and intensity of the Devonian field reported here, although it is not yet clear whether erratic directional behaviour should always correlate with weak palaeointensities and vice versa. Various mechanisms for producing transitions between the strong-field dipole dominated state and these states, linked to surface manifestations of mantle convection, have been proposed and include superplume growth and collapse (Amit & Olson 2015) and episodes of true polar wander (Biggin et al. 2012). Yet another mechanism would be inner core nucleation (Driscoll 2016) but this would require its delay until after 400 Ma in order that the weak-field multipolar dynamo state can still be observed in the Devonian.
If the field was in a multipolar state for long periods, then the resulting diminished magnetosphere would have increased the impact of solar particles on the Earth's ionosphere and atmosphere (Siscoe & Chen 1975;Siscoe & Crooker 1976;Starchenko & Shcherbakov 1991;Vogt et al. 2007). The characteristic dayside size R m of the magnetosphere can be estimated by balancing the solar wind pressure on the sunward side of the boundary and the magnetic pressure on the earthward side (Martyn 1951;Siscoe & Chen 1975) yielding the scaling relation R m = const × B s n where B s is the intensity of the field on the Earth surface and the exponent n = 3 for dipolar configuration, n = 4 for the quadrupolar configuration and so on. Assuming a dipole-dominated field, a tenfold decrease of the intensity will diminish the dayside size by about a factor of two from its present-day value R m ≈ 10 R E where R E is the Earth's radius. This would drastically affect the geomagnetic shielding efficiency and bring substantial changes to the Earth's ionosphere and atmosphere. In particular, Starchenko & Shcherbakov (1991) and Vogt et al. (2007) demonstrated that a dangerous situation involving penetration of solar protons of several tens of MeV into Earth's atmosphere arises when patches of high magnetic inclination are found at low latitudes opening new wide particle-entry regions. The possible impacts of enhanced solar wind penetration on the climate and biota of the Devonian warrant further investigations and testing against observational records.
For the presentation of a general view on the problem of possible complex field geometry, Fig. 15 displays PINT data not only for the Devonian but also for the neighbouring Silurian and Carboniferous geological periods. Remarkably, reported data for the early  The VADMs for the Minusa Basin results were calculated using the inclination derived from the N and S components (41 • ) and those for the Kola results were calculated using the average inclination for the D component (16 • ) provided by Veselovskiy et al. (2013). For explanation of the Q PI criteria, readers are directed to the main text and .
Carboniferous also reveal low and very low intensities. There are only two data available for the Silurian and both of them show low B anc , These observations hint that the epoch of unstable geomagnetic field generation regime could be extended for the time interval somewhat wider than just the Devonian. Unfortunately, poor statistics and the fact that all the data published before do not satisfy the modern reliability criteria makes such the suggestions speculative. In fact, our analysis have clearly shown that the Palaeozoic palaeomagnetic data are much too far sparse and a number of new high quality palaeointensity and palaeodirectional determination are required to make definite conclusions on the characteristics of geomagnetic field for these extremely interesting geological periods.

C O N C L U S I O N S
We have reported evidence that the time-averaged geomagnetic field between ca. 375 and 408 Ma was considerably weaker than that of the last few millions of years, or indeed of that encountered at any point in the last 400 Ma except perhaps the mid-Jurassic. Our estimates of the virtual dipole moment are drawn from 12 sites from two regions on distinct cratons argued to be separated by more than 15 • in palaeolatitude and likely with a substantial body of water in between them . Within-site and between-method variability suggests that our palaeointensity results are imprecise on a relative scale. Nevertheless, they are near-uniformly low in value (producing VDMs that are, except in two cases, all below 20 ZAm 2 ), display a high degree of internal consistency on an absolute scale (site mean standard deviations generally less than 5 µT), and have high associated Q PI values.
Together with the peculiar directional behaviour of the Devonian field recognized in earlier studies, these findings suggest that the time-averaged field may have been in a state of far reduced dipole dominance and perhaps even in a persistently multipolar state. The duration, extent, and precise nature of this state await further investigations, some of which are already underway. The implications for the geodynamo, core-mantle interaction, the Devonian environment and evolutionary pressures on biota are currently unknown, but are likely to be significant.

A C K N O W L E D G E M E N T S
The work was supported by Russian Fund of Basic Research (grant nos 15-35-20599, 15-55-10055KO-a, 16-05-446-a), Ministry of Education and Science of the Russian Federation (grant 14.Z50.31.0017). AJB acknowledges funding from the Natural Environment Research Council (NE/P00170X/1) and The Royal Society ( IE150283). LMAH acknowledges funding from a NERC DTP Studentship (1511981). The authors are grateful to Rob van der Voo and Hidetoshi Shibuya for the valuable comments and suggestions made to this paper. We thank Vladimir Pavlov for fruitful discussion and Vladimir Tselmovich for the help in producing the microprobe studies.

S U P P O RT I N G I N F O R M AT I O N
Supplementary data are available at GJIRAS online. Table S1. Palaeointensity results at the specimen level. Shaded rows indicate rejected samples with the failing criteria underlined. Methods used were Microwave IZZI protocol with pTRM checks (MW), Microwave Coe protocol with pTRM and tail checks (MWC), thermal Thellier Coe protocol with pTRM checks and occasionally tail checks (TT/TT(2c)/TT(BL) -note 2c refers to where the sample magnetization was only measured in two components, BL refers to baseline where measurements were made above room temperature). For MW results, the applied power in Watts and the duration of the application in seconds is given in place of sample peak temperature. All parameters are described in the standard palaeointensity descriptions (SPD; Paterson et al. 2014).
Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.