Constraining the eruption history of Rangitoto volcano, New Zealand, using palaeomagnetic data

is situated within the


Introduction
Small basaltic volcanic edifices located on continental crust commonly occur in clusters or fields.The individual volcanic vents associated with the edifices are referred to as monogenetic as they are thought to erupt only once, and the activity is short-lived (weeks to months) (Walker 1993).The eruption recurrence rates at some volcanic fields can be thousands of years so it is difficult to assess when or where the next eruption will occur and how long it might last.
Rangitoto Island Volcano (hereon referred to as Rangitoto) is a shield volcano located in the Waitemata Harbour of Auckland city, New Zealand (Fig. 1).The volcano is part of the basaltic and monogenetic Auckland Volcanic Field (AVF), which has been active over the past ca.250,000 years (Kermode 1992) with ~50 eruption centers in the form of maars, scoria cones and tuff rings (Hopkins et al., 2021).Rangitoto is the youngest and, by far the largest, volcano in the AVF, suggesting that the AVF may have entered a new and different phase of activity (Linnell et al., 2016).Given the proximity of Rangitoto to the city of Auckland, its eruption history is important for assessing future eruption style and duration scenarios for hazard and risk planning.
The age of the Rangitoto lava shield is uncertain because the lava flows have not been directly dated.One hypothesis is that the volcano, including the shield and summit scoria cones, formed in one short episode (i.e.-monogenetic volcanism).Lindsay et al. (2011) suggest an age of around 500-550 cal yr BP for the entire volcano based on radiocarbon ages on organic material interbedded with tephra on adjacent Motutapu Island (Fig. 1C).However, it has not been demonstrated that the tephra was synchronous with the lava effusion and might relate only to the volcanism that produced the overlying scoria cones.Furthermore, Shane et al. (2013) highlighted that some of the radiocarbon ages are stratigraphically inverted and include age determinations of up to 890 cal yr BP from those above the tephra layers.Based on crypto-tephra in lake sediment, Shane et al. (2013) identified possible Rangitoto activity back to 1498 ± 140 cal yr BP, representing a total lifespan of approximately 1000 years.Linnell et al. (2016) proposed early volcanism may have started at 6 ka based on a single radiocarbon dated lava flow interbedded with sediments beneath the shield edifice.Alternatively, this flow could represent an earlier, buried volcano.On this basis, Linnell et al. (2016) proposed a timespan of ~650-550 cal yr BP for the shield building phase based on the youngest radiocarbon age from the sediments underlying the lava shield and using the age of the 550 cal yr BP from Lindsay et al. (2011) representing the start of the pyroclastic activity (scoria cones) that overlie the shield of Rangitoto.
In a previous palaeomagnetic investigation of the Rangitoto lava flows, Robertson (1986) measured palaeodirections from eleven sites on the flanks of the volcano.If the lavas erupted over a very short timescale, it would be expected that the magnetic directions recorded by the lavas would show little variation (due to the expected slow rate of change of the geomagnetic field).The average declination for the sample sites varied between ~354 • and ~13 • and the inclinations varied between approximately − 56 • and − 71 • .Assuming a constant rate of change in declination, which we now recognise to be a reasonable approximation for the past ~700 years in New Zealand (Turner et al., 2015a), Robertson (1986) argued for eruption activity over 1000 years.Overall, the scatter between the site mean directions, even considering the rather large uncertainties, is incompatible with a short duration for the construction of the shield volcano, unless the lava flow carapaces partly congealed prior to final emplacement.The absence of an unambiguous chronology for the Rangitoto volcano shield-building phase, and in particular as direct radiocarbon ages on eruption deposits have evidence for reworking (see Shane et al.,

2013
), drove the present investigation into the potential of the palaeomagnetic record to provide new constraints on the eruption history using lava flows sampled by deep drilling in the flank of Rangitoto volcano.The radiocarbon ages used to constrain the age of the drill core have also been re-evaluated to incorporate the uncertainties from the marine reservoir effect (Petchey and Schmid 2020), which had previously not been considered.Magnetic mineralogy experiments were also undertaken for strategically selected samples from the lava flows to verify that the lavas are suitable recorders of the geomagnetic field.Finally, the minimum duration of the shield-building phase required by the palaeomagnetic data has been evaluated using a statistical modelling approach.

Geological setting
The AVF (Fig. 1) is a Quaternary basaltic volcanic field in a continental intraplate setting that has a crustal thickness of approximately ~25-30 km (Horspool et al., 2006).The field covers an area of 360 km 2 and consists of around 50 basaltic cinder cones, maars and associated lava flows (e.g.-Kereszturi et al., 2014).Rangitoto volcano (36.79 • S, 174.86 • E) is a near symmetrical volcano 6 km in diameter, a maximum elevation of 260 m and a volume of approximately 0.7 km 3 (Kereszturi et al., 2014).The volcano is significantly larger than the others in the AVF and is responsible for just under half of the total erupted magma in the field (Linnell et al., 2016).The summit of Rangitoto consists of a central scoria cone flanked by remnant scoria mounds and ridges from older cones to the north and the south.Geochemical analysis of the surface lava flows has shown that the central and southern cones erupted subalkalic basalts whilst the north cone produced alkalic basalt (McGee et al. 2011;Needham et al. 2011;McGee et al. 2013).

Fig. 2.
Lithological log of Rangitoto drill core, modified from Linnell et al. (2016), showing the different identified lava flows and the geochemical zones.The ages shown are the uncalibrated radiocarbon ages, also from Linnell et al. (2016).The open circles show the location of the core samples.

Methods
The samples for the palaeomagnetic experiments described here were obtained from the Rangitoto volcano core recovered and first described in Linnell et al. (2016).The ~150 m vertically-drilled core was taken from the western flank of Rangitoto at metric grid north 5927446, east 1765100 (map NZTopo50-BA32, edition 1.02; Land Information, New Zealand, 2012) at approximately 120 m a.s.l.(Fig. 1B).
The top 128 m of the core is composed of at least 49 distinct massive basaltic lava flow units, separated by autobreccias, which are interpreted to represent the main shield building phase (Fig. 2).These flows are underlain by (i) ~15 m of marine sediments interbedded with thick basaltic volcanic or pyroclastic deposits, (ii) an older (~60 cm thick) lava flow (dated to ~ 6000 cal years BP) and (iii) Miocene deep marine sandstone at the base of the core.Linnell et al. (2016) produced a geochemical stratigraphy for the volcano based on the succession of flows in the drill core with four distinguished subalkalic basalt magma groups: group 1 (137-95 m depth), group 2 (94-69 m depth), group 3 (68-26 m depth) and group 4 (26-0 m depth), which are also shown in Fig. 2. Thirteen radiocarbon dates were taken from sediments below the lava flows from the shield building part of the core.As can be seen in Fig. 2, several of these ages are inverted suggesting substantial reworking of the sediments.
The core samples were obtained from one continuous near-vertical hole.The recovered intact core sections were not oriented with respect to geographic north, which means that it was only possible to retrieve palaeoinclinations, i.e.-the dip angle of the magnetisation vector.As all the palaeomagnetic samples were taken from the same core, they are not independently orientated and increasing the number of samples will therefore not necessarily reduce the total error (not considering orientation errors).
Preliminary demagnetisation experiments showed that the upper sections of the core (above 115 m) carried more stable (higher coercivity) magnetisations.Due to availability of samples, lava flows 1-28 (between 0 and 70 m) were sampled more heavily, but all depths and lava flows that could be sampled are represented to some degree.The subsamples consist of standard size (2.5 cm diameter) palaeomagnetic cores, drilled perpendicular to the length of the drill core.The short cores were subsequently cut into smaller (approximately 1 cm 3 ) subsamples for the different palaeomagnetic and mineral magnetic experiments.The subsamples were taken from the centre of the drill core in order to minimise any effects of a drilling induced remanence.
Hysteresis loops and temperature dependence of magnetisation were measured for 9 samples (spaced between 8.20 m and 127.95 m) using a Variable Field Translation Balance (VFTB) at the Geomagnetism Laboratory, University of Liverpool, UK.Sample preparation involved grinding approximately 150 mg of each sample to a fine powder which was then transferred to an appropriate sample holder.The heating experiments were conducted in air between 40 • C and 700 • C. The VFTB results were analysed using the software RockMagAnalyser1.0 (Leonhardt 2006).The temperature dependence of magnetic susceptibility was measured for 24 samples (from a range of depths between 7.80 m (lava flow 1) and 127.95 m (lava flow 50) on an MFK1-FA Kappabridge fitted with a CS4 attachment at the Palaeomagnetic Laboratory, Lund University, Sweden.Sample preparation involves grinding approximately 0.25 cm 3 of the sample to a fine powder using a mortar and pestle.The experiments were completed in an argon atmosphere between 40 • C and 700 • C. To further investigate the alteration at mid-to high-temperatures selected samples were also subjected to temperature cycling experiments in which the samples were heated to progressively higher peak temperatures.
The magnetic remanence directions were obtained using a 2G Enterprises cryogenic magnetometer (model 760) equipped with an automated alternating field (AF) demagnetisation system at the Palaeomagnetic Laboratory, Lund University.The natural remanent magnetisation (NRM) of 157 samples were measured before and after AF demagnetisation in 5 mT steps up to 40 mT, then in 10 mT steps until a peak field of 80 mT was applied.Replicate palaeodirections were also measured for 37 of the palaeomagnetic subsamples to check for consistency.The palaeodirections were calculated using principal component analysis (Kirschvink 1980) based on the NRM data between 15 and 80 mT.
Palaeointensity experiments were undertaken at Lund University using the 2G 760 magnetometer and a Magnetic Moments Thermal Demagnetiser (MMTD80) and at the University of Liverpool using a microwave system (MWS) consisting of a 14 GHz microwave system combined with a low temperature Tristan SQUID magnetometer.Microwave and thermal methods for determining palaeointensities have been demonstrated to provide comparable results (e.g.-Hill et al., 2002).All the palaeointensity experiments were run using the IZZI + protocol (Yu et al., 2004) and a laboratory field of 50 μT.54 samples (distributed throughout lava flow 1 to lava flow 49) were analysed with the MWS and 8 samples with conventional Thellier (samples every 10 m downcore between 30 and 100 m, which corresponds to flows 14-40).Pseudo-Thellier experiments, initially designed for sediments (Tauxe et al., 1995) were also undertaken on 37 samples spaced along the whole length of the core which contains lava.Using the 2G 760 magnetometer the selected samples were first stepwise AF demagnetised using the same steps as detailed above and then given an anhysteretic remanent magnetisation (ARM) with successively higher peak alternating fields and using a direct current bias field of 50 μT.

1: Magnetic mineralogy
To assess the suitability of the Rangitoto lavas for palaeomagnetic experiments, a suite of magnetic mineralogy experiments was undertaken.Overall, the mineralogy results are somewhat typical for basaltic lava flows (e.g.-Petersen 1976) indicating the presence of (titano) magnetite with varying proportions of titanium with a range of oxidation having occurred (Fig. 3).Approximately 115 m downcore, the mineralogy changes significantly with lavas below this depth, interpreted to be dominated by high titanium titanomagnetite.We note, that this transition does not appear to be related to any major changes in bulk geochemistry of the lava flows (Linnell et al., 2016) (see also Fig. 2).
The results from the hysteresis measurements have been summarised in a Day plot (Day et al., 1977) (Fig. 3A), which indicates that the nine samples measured (covering all four geochemical zones) are dominated by bulk pseudo-single domain grain sizes as is typical for lava flows.The samples generally saturate around applied fields in the range of 300-400 mT.Hysteresis loops from samples above 115 m (Fig. 3B) are characterised by relatively high coercivities (B C > 20 mT), whilst hysteresis loops from samples below 115 m (Fig. 3C) have a narrower profile, indicating lower coercivities (B C < 10 mT).
Temperature dependence of magnetisation for samples above 115 m show a characteristic drop in intensity around 580 • C upon heating (Fig. 3D), indicating that magnetite is a main magnetic carrier, and slightly lower intensities upon cooling suggesting minor degrees of alteration.For samples below 115 m, the intensity drops rapidly around 100 • C upon heating (Fig. 3E).We interpret this as high-titanium titanomagnetite, which oxidises into a much stronger magnetic phase (magnetite) when heated in air (Fig. 3E), but remains stable when heated in an argon atmosphere (Fig. 3F).Based on this interpretation, we anticipate that the lavas from below 115 m could provide suitable recordings of past geomagnetic field directions, however, it is expected that there may be limited success in the palaeointensity experiments for those samples susceptible to alteration on heating in air.
Fig. 3G -I) shows a selection of results from the CS4 magnetic mineralogy experiments for samples above 115 m.The characteristic Hopkinson peak for SD magnetite is present, supporting the interpretation that magnetite is a main magnetic carrier for the majority of these  (Linnell et al., 2016).Purple = geochemical group 1 (bottom of core), yellow = group 2, orange = group 3 and blue = group 4 (top of core).Day Plots (Day et al., 1977) plot the ratio of saturation remanent magnetisation to saturation magnetisation (Mrs/Ms) versus the ratio of coercivity of remanence to coercivity (Bcr/Bc).SD = single domain, PSD = pseudo-single domain and MD = multi domain.Theoretical SD/MD mixing curves (dashed lines) from Dunlop (2002aDunlop ( , 2002b) ) and 460 • C (41.50 m) and (L) between 500 • C and 590 • C (3.73 m).samples.Alongside magnetite, a number of other magnetic phases were identified and thus the temperature cycling experiments were tailored to check for potential alteration associated with these different magnetic phases.Fig. 3J shows a cycle investigating a low temperature magnetic phase (100-150 • C) in this case from a sample at 3.73 m (top of flow 1).There is clear evidence for reversible behaviour between 90 • C and 170 • C, suggesting that the peak and associated drop in magnetic susceptibility around 100-150 • C is a Curie/Neél temperature.Fig. 3K investigates the magnetic phase between 320 • C to 460 • C.There is evidence for some degree of irreversibility, leading to slightly higher magnetic susceptibilities with heating to higher temperatures.However, the drop in magnetic susceptibility remains stable at around 350 • C indicating that this is likely also a Curie/Neél temperature (likely titanomagnetite).Fig. 3L investigates potential alteration occurring between 500 • C and 590 • C. In this case, the drop in magnetic susceptibility between ~500 • C and ~550 • C is clearly associated with irreversible behaviour indicating thermochemical alteration of this magnetic phase.This is a more extreme example of alteration in the upper 115 m of the core, the majority of samples show significantly less evidence of alteration than shown here.In summary, the temperature cycling experiments (Fig. 3J-L) show that the magnetic mineralogy is stable with heating up to 500 • C (and reversible in argon) after which irreversible behaviour is observed in at least two of the cycles which could be an indication for the presence of maghaemite (partially) inverting to haematite (Dunlop and Özdemir 1997).If maghaemite is present as a result of oxidation, then those samples are not suitable for palaeointensity experiments which require a thermal remanence.

Palaeomagnetic inclinations
With the exception of a few samples from the middle of flow 1, the NRM is carried by a single component magnetisation that demagnetises (AF, thermal, microwave) to the origin after removing a small overprint (Fig. 4A-H).NRM intensities varies from approximately 1 to 90 Am − 1 with a median value of 12 Am − 1 , all calculated assuming a volume of 1 cm 3 .The lack of a significant drilling overprint is likely due to the sampling strategy as samples were taken from the middle of the core approximately 4-5 cm from the core edge.Based on the uncomplicated demagnetisation behaviour, we focus solely on directional results obtained using AF demagnetisation (Fig. 5).
Palaeomagnetic inclinations vary from − 75 • to − 40 • (Fig. 5A), encompassing the inclination of − 56.2 • calculated at this location for a geocentric axial dipole (GAD) field and the present-day inclination of − 62.7 • .Maximum angular deviations (MAD) are predominantly below 1 • (Fig. 5B) and angular deviations from the origin (ADdefined as the angle between a free fit and anchored PCA fit) are mostly below 2 • (Fig. 5C), indicating a stable characteristic remanent magnetisation.Downcore variations in the median destructive field of NRM (MDF NRM ) (Fig. 5D) are consistent with the hysteresis results (Fig. 3), with generally high values (20-60 mT) above 115 m and low values (15-20 mT) below this depth.Due to the high coercivities of some samples (MDF NRM > 70 mT), the peak AF field (80 mT) was sometimes not sufficient to fully demagnetise the NRM, but the low AD values show that the isolated component is still heading to the origin.Notably shallow inclinations are recorded around 5-8 m core depth (lava flow 1) and between 36 and 45 m depth (lava flows 16-20).In the case of lava flow 1, we note that the shallower inclinations are concentrated in the centre of the flow and coincide with an increase in MAD values (from under 1 • -2.5 • ) and increased AD values (up to 7 • ), perhaps indicating rheological deformation of the flow during cooling.The opposite is true for flows 16-20 which represent some of the best defined directions (Fig. 5B-C).The calculated mean (per flow) does not change substantially if the inclination from the centre of flow 1 (5-8.5 m depth) are included (− 54.8 • ± 1 • ) or excluded (− 57.1 • ± 1.1 • ) when compared to the scale of variations seen throughout the core and we therefore opted to not exclude any of the results.

palaeointensity results
Palaeointensity experiments were carried out using a combination of the pseudo-Thellier approach (Fig. 4I), conventional Thellier (Fig. 4J) and MWS (Fig. 4K-L).Absolute palaeointensity experiments were carried out on 62 samples (8 conventional Thellier and 54 MWS).Using the SELCRIT2 criteria (Biggin et al., 2007;Paterson et al., 2014) we obtained a total of 33 successful absolute palaeointensity estimates.The primary reason for rejecting palaeointensity estimates was due to failed pTRM checks indicating thermochemical alteration, particularly at high temperatures/high microwave energy input as was expected given the magnetic mineral investigation (Fig. 3).The palaeointensity experiments that were conducted on the MWS-which was designed to minimise thermochemical alteration (Hill and Shaw 1999;Suttie et al., 2010), had a higher success rate (61%) when compared to the conventional palaeointensity experiments (12.5%, or only 1 of 8) and an overall success rate of 53%.An additional factor may be the high unblocking temperatures of the samples (Fig. 4J) making the experiments sensitive to the reproducibility of the oven temperature.
The accepted absolute palaeointensity estimates show consistent downcore variations with values around 40-50 μT in the top of the core, a minimum, potentially as low as ~25 μT around 20 m depth and values around 50-55 μT between 40 m and 60 m (Fig. 5E).Below 60 m, there are too few data (three scattered estimates between 35 and 55 μT) to determine any trends.The intensity variations in the top 40 m of the core appear to be somewhat synchronous with the changes seen in the inclination values, whereby the low intensities around 20 m coincide with the steeper inclinations.This supports the hypothesis that these variations are geomagnetic in origin, although we acknowledge that steeper inclinations are typically associated with stronger intensities in a dipole dominated field.The pseudo-Thellier data (Fig. 5F) estimates can be calibrated to absolute values (F ABS ) using a linear relationship (given below) assuming that there are no substantial variations in magnetic grain sizes (or mineralogy and concentrations) (Tauxe et al., 1995).A common approach is to use the median acquisition field of ARM (MAF ARM ) (Fig. 5D) as a proxy check for variations in grain size and to reject samples that fall outside a certain range, e.g.-between 23 and 63 mT (de Groot et al., 2013).Based on these criteria, 8/47 samples from the pseudo-Thellier experiments were rejected.The remaining MAF ARM values do not show substantial variations (Fig. 5D) suggesting that the observed relative palaeointensity changes can indeed be calibrated to absolute values.Here, we use the linear function F ABS = 7.718 × | pseudo-Thellier slope| + 14.600 (de Groot et al., 2016), but note that this function will likely lead to a slight overestimation of the intensities since it assumes a DC bias field of 40 μT (whilst these experiments use 50 μT).Overall, there is good agreement between the absolute palaeointensity data and the pseudo-Thellier results (Fig. 5E) with relatively high palaeointensities between 40 and 60 m in the core and low intensities around 20 m depth in the core.However, we also note that the results diverge for the top of the core so that it is not possible to rule out that part of the variations observed in the pseudo-Thellier data (and to a lesser degree the absolute intensity data) are related to changes in the mineralogy associated with the different geochemical zones.

Reinterpreting the radiocarbon ages
Linnell et al. ( 2016) placed the main shield-building phase between 650 and 550 cal years BP.The upper bound (approximately 550 cal years BP) was determined based on the assumption that the radiocarbon dated tephra layers from nearby Motutapu swamp sediment cores are associated with pyroclastic activity from the summit scoria cones.Shane et al. (2013) note that some radiocarbon ages associated with the tephra layers are stratigraphically inverted and that reworking may have occurred.This places additional uncertainty on the upper age boundary.The lower bound (approximately 650 cal years BP) is based on the youngest radiocarbon age (1092 ± 26 14 C years) from foraminifera extracted from the marine sediment below 128 m in the core.However, some radiocarbon ages from various macrofossils found beneath this horizon are also inverted which indicates reworking of sediments has occurred (Linnell et al., 2016).The calibrated ages for the youngest radiocarbon date (~650 cal years BP) from the sediment beneath the lava was obtained originally using the Marine13 calibration curve (Reimer et al., 2013) assuming an average marine reservoir effect (local reservoir correction ΔR = 0).For the sake of continuity, all calibrations have been done using SHCal13 (Hogg et al., 2013), to match the analysis undertaken on the core in Linnell et al. (2016).
Due to the proximity of Rangitoto to the Auckland coastline (less than 10 km from the coast) it is reasonable to expect some freshwater input (with atmospheric 14 C concentrations), which suggests that the local marine reservoir age might be lower than average marine reservoir age representative of the oceanic 14 C concentrations (e.g.-Lougheed et al., 2013).Such local deviations from the average marine reservoir age (and associated uncertainties) are typically accounted for by introducing a local reservoir correction, ΔR.The true calendar age of the youngest radiocarbon date obtained from the sediments below 128 m of the core likely falls between two end members: (i) assuming a marine signal with ΔR = 0 14 C years yielding a minimum age of ~650 cal years BP as in Linnell et al. ( 2016) and (ii) assuming a purely atmospheric signal and using the SHCal13 calibration curve (Hogg et al., 2013), yielding a maximum probability age of ~950 cal years BP.The latter case is approximately the equivalent of using the Marine13 calibration curve and a ΔR = − 325 14 C years (Supplementary Fig. S1).Reconstructions of variations in the marine reservoir effect in the Hauraki Gulf region based on independently dated samples (Petchey and Schmid 2020) shows large variations in ΔR over time, particularly around the period 550 and 600 cal years BP, with the calculated mean ΔR reaching as low as − 172 ± 92 14 C years.Based on a detailed investigation of the compiled ΔR values within the Hauraki gulf embayment, spanning the time period 300-750 cal years BP (Fig. 6), we propose a more appropriate local mean ΔR value of − 91 ± 139 14 C years.Using the Marine13 calibration curve and this updated ΔR value to calibrate the youngest radiocarbon date (from below the lava flows) we get ages with a 95% range of 523-1028 cal years BP.This result is consistent with the two endmembers proposed above and suggests an older age for the lower bound of the shield-building phase, implying possible total durations of up to 400-500 years.
Acknowledging the risk that any radiocarbon determination could suffer from contamination (i.e., providing incorrect ages), we construct four potential scenarios (I -IV) for the total duration of the shieldbuilding phase (T) which we will refer to in the continued discussion.For scenario I and II we assume a fixed upper bound of 550 years BP (from the dated scoria cones).The lower bounds are based on the calibrated ages of each of the two youngest radiocarbon dates from the core (1092 ± 26 and 1429 ± 20 14 C years BP).Using the two end-member argument formulated above, we define a maximum duration T MAX = 505 years for scenario I and T MAX = 801 years for scenario II based on the 99.9th percentile of the calibrated ages using SHCal13 (Supplementary Fig. S1), i.e., assuming a purely atmospheric signal.For scenario III, we relax the upper age constraint to 350 cal years BP (based on palaeomagnetic model-data comparisons, see section 5.2) and set an arbitrary large T MAX = 2000 years which in effect allows the age determination to be influenced solely by the palaeomagnetic observations.Finally, we add a fourth scenario (IV) with the same relaxed upper age constraint of 350 cal years BP and the lower bound defined as in scenario I. Table 1 shows an overview of the four scenarios.

Geomagnetic model-data comparison
The new palaeomagnetic results for Rangitoto (Fig. 7A and B) can be compared to independent archaeomagnetic data (burnt archaeological artefacts and volcanic rocks) from New Zealand, sourced from GEO-MAGIA50.v3 (Brown et al., 2015), the palaeosecular variation (PSV) record from Lake Mavora, New Zealand (Turner et al., 2015b), the NZPSV1k master curve (Turner et al. 2015a(Turner et al. , 2020) ) and the predictions of two geomagnetic field models: pfm9k.2(Nilsson et al., 2022) and COV-LAKE (Hellio and Gillet 2018) (Fig. 7C and D).For the comparison we calculated the average inclination per flow (where three or more inclination values were available) and used the Student's t-distribution (the one-dimensional equivalent of Fisher's cone of confidence) to determine uncertainties (Fig. 7A).For the intensity data, there are generally not enough datapoints to determine the uncertainties and we therefore only show averages for the absolute and pseudo-Thellier estimates (Fig. 7B).
The global geomagnetic field models pfm9k.2(Nilsson et al., 2022) and COV-LAKE (Hellio and Gillet 2018) are constrained by a combination of archaeomagnetic data and palaeomagnetic data from lake/marine sediments.Both models use a similar statistical approach to incorporate prior information from satellite era models to constrain the potential variability of the geomagnetic field and indirectly the model uncertainties (represented by different model realizations in Fig. 7C and  D).The differences in model predictions are mostly due to the selection of data used to constrain the models, and to a lesser degree due to the modelling strategies (primarily the treatment of age uncertainties and sedimentary data).As shown in Fig. 7C and D, there are only a limited amount of archaeomagnetic data available from the region and the COV-LAKE model prediction is therefore largely influenced by the sedimentary data-primarily the PSV record from Lake Mavora (which is also shown in Fig. 7C and D).The NZPSV1k master curve contains data from Lake Mavora and Lake Ponui and also includes the gufm1 historical field model for the past 400 years (Jackson et al., 2000), although we note that inclination observations from the region are only available from the late 18th century.Due to the gradual process by which sediments acquire a magnetisation, models that include sedimentary data tend to be complicated by unaccounted for temporal smoothing and time lags in the data (Mellström et al., 2015;Nilsson et al., 2018;Nilsson and Suttie 2021).These undesirable effects are what pfm9k.2 was designed to avoid, but due to the computational cost associated with each new sediment record, Lake Mavora and other records from New Zealand were not included in this version of the model.This omission is partially accounted for by the larger dispersion of the model predictions.Nevertheless, the global geomagnetic field models have the advantage of appropriately incorporating information from more distant data sources (e.g.-from Australia) that provide constraints on the large-scale field variations.
Large and statistically significant changes in inclination between flows are clearly visible in the Rangitoto data (Fig. 7A), in some cases by as much as 10 • -15 • (seen between flow 15 and 16).Whilst there are large differences between the reference model predictions (and datasets) (Fig. 7C), it is clear that such large swings in inclination are normally associated with geomagnetic field variations on centennial (or even millennial) timescales.The closest analogue (both in terms of variability and timing) to the recorded variations, in particular the shallow inclinations between flow 16 and 20, is the field prediction by the COV-LAKE model over the period 1500 to 500 years BP.This is consistent with scenario III, and potentially also scenario II, but definitely inconsistent with scenario I.The gradual shallowing in inclination between flow 10 and 1 also looks strikingly similar to the inclination variations predicted by both global models between 550 and 350 years BP, but less apparent in the Lake Mavora record and absent in the NZPSV1k master curve.The similarity between the model and data is perhaps indicating that the proposed upper bound of the shield building phase (~550 years BP) dated by the appearance of summit scoria cones may need to be revised.As observed by Shane et al. (2013), some of the radiocarbon ages from the tephra layers identified in the Motutapu swamp are not consistent with stratigraphic order, suggesting that the timing of the pyroclastic activity is likely not well-constrained.This apparent uncertainty in the upper age of the shield building phase is the motivation for introducing scenario IV (see Table 1).
We note an overall trend of stable intensities around 50-60 μT followed by a period of slightly weaker intensities (30-40 μT) between flows 7 to 15, with a return to higher intensities in the top of the core (flows 1-6) (Fig. 7B).These variations are generally consistent with the COV-LAKE model predictions between 1500 and 500 years BP, and perhaps even more so with Lake Mavora PSV record over the same time period (Fig. 7D).We note, however, that a satisfactory theoretical foundation for palaeointensity estimation of sediments is still lacking (Roberts et al., 2013) and the palaeointensities recorded by Lake Mavora (that are also used to constrain COV-LAKE) should therefore be treated with caution.
Overall, we note that there are large disagreements between the different reference curves and the reference data (Fig. 7C-D), which limits the possibilities to constrain the Rangitoto eruption history using palaeomagnetic data.It is, however, clear from the model-data comparison that the reconstructed variations in both inclination and intensity from Rangitoto are more consistent with a relatively long duration for the shield building phase (e.g.-scenario III) rather than a short duration (scenario I).This implies that the observed palaeomagnetic variations (both inclinations and intensities) are inconsistent with the radiocarbon data and that additional, and currently unaccounted for, sources of uncertainty may need to be invoked to explain some of these observed variations.Processes that could produce additional sources of uncertainty can be broadly separated into two categories; those that affect the physical orientation of the samples and those that affect the chemical properties of the lavas.
The first category includes processes such as tectonic tilting, rheological deformation and/or the angle at which the coring drill penetrated the lavas, all of which could affect inclination values.Posteruption tilting of the lava flows could potentially explain the shallow directions recorded in flows 16-20 (Speranza et al., 2012) but we note that there is no evidence of tectonic activity in the Auckland area over the last few thousand years (Kermode 1992).Rheological deformation of the lava flows, i.e., when the outer carapace become solid while the interiors remain molten, could also affect the palaeomagnetic directions.Any such deformation (or other movement) occurring when the sampled part of the lava flow has cooled below the TRM blocking temperatures could result in single component magnetisations (as is the case for the Rangitoto samples), with systematically too shallow or too steep inclinations ( de Groot et al., 2016).Normally, such effects would be detected by sampling different outcrops of the same lava flow, which is not possible with the drill core sampled in this study.If the core drilling and recovery process allowed for the possibility that lava flows were sampled at an angle other than vertical, it would likely be a source of some systematic uncertainty.However, liaison with the drill team have confirmed that the drill hole was near vertical (otherwise the rigid drill rods would have become stuck), which suggests that coring angle is not a large source of uncertainty in this study.This interpretation is supported by the overall good agreement between the inclinations obtained in the present study and the palaeodirections obtained from 11 flows by Robertson (1986), with a site mean inclination 61.7 The original TRM could also have been replaced by a hard IRM obtained from an event like a lightning strike (Di Chiara et al., 2012), which would cause errant palaeodirectional and palaeointensity values.Geochemical alteration (due to weathering processes etc.), will mainly affect the obtained palaeointensity values.Maghaematisation does not generally affect the direction of the remanent magnetisation (Dunlop and Özdemir 1997).The rock magnetic analysis has indicated the possible presence of maghaemite in some samples, which suggests a variable degree of oxidation through weathering.This suspected presence of maghaemite plus the variable degree of thermochemical alteration observed in the thermomagnetic curves (Fig. 3) suggests that the palaeointensity data should be treated with caution.While we cannot rule out other factors that may influence the directional data (e.g.magnetic refraction, flow anisotropy) apart from the middle of flow 1, there are no obvious indicators that the data has been affected in some way.On the contrary, we note that many of the observed variations in inclination and intensity are coeval, thus supporting that these variations represent changes in the geomagnetic field.

Modelling
Due to the large uncertainties in the model reference curves (large disagreements between different reference datasets), it is difficult to use the palaeomagnetic data to directly date the Rangitoto shield building phase (Fig. 7).Instead, we propose a Bayesian statistical model to constrain the duration of the shield building phase using prior information on the expected rates of change of the geomagnetic field to estimate the minimum time required to explain the variations seen in the palaeomagnetic data.Here we provide a brief outline of the model, with a detailed description given in the supplementary material.Similar to Nilsson et al. (2022) and Hellio and Gillet (2018), we assume that the geomagnetic field results from a stationary stochastic process and not all flows are equal in size and some flows are more heavily represented in the material received for sampling.(B) Average absolute (black) and calibrated pseudo-Thellier (grey) palaeointensities per lava flow.(C) and D) independent inclination and intensity data from New Zealand, relocated to Rangitoto assuming a geocentric dipolar field: archaeomagnetic data (black) from Geomagia50.v3 with 2-sigma error bars, the Lake Mavora PSV record with 2-sigma uncertainty envelope (red) and NZPSV1k with 2-sigma uncertainty envelope (cyan).The archaeomagnetic data from Geomagia50.v3 without independent age control or with age uncertainties in excess of 2000 years were removed.The data from Turner et al. (2020) were re-plotted using the independent ages provided in the original publication.The relative palaeointensity data from Lake Mavora were converted to absolute values using the scaling factor proposed by Turner et al. (2015b) prior to being relocated.Also shown are the predicted inclinations and intensities of two geomagnetic field models: pfm9k.2(blue) and COV-LAKE (green), with the mean model prediction shown as bold lines and 50 possible realizations as thin semi-transparent lines.The four scenarios (I -IV) for the duration of the shield building phase are represented by arrows above (C) and (D).
determine its statistical properties from satellite era models (Nilsson and Suttie 2021).We note that this prior description of the field does not preclude the possibility of variations more rapid than what has been observed in the satellite-era if the data requires it.The model has been designed to estimate geomagnetic field variations in inclination only, using the inclination data from the N = 24 lava flows that can be considered to have well defined uncertainties (Fig. 7A).The total duration of the shield building phase is defined as T = ∑ N i=1 Δt i , where Δt 1 , …, Δt N− 1 represents the amount of time that has passed between each of the N sampled lava flows.The likelihood function comes down to evaluating the probability of the observed variations in inclination, given a certain duration of the shield building phase, T, which determines the range of "allowed" variability in the modelled inclination time-series.To sample from the posterior distribution, we use a Markov Chain Monte Carlo (MCMC) approach, specifically the No-U-Turn sampler (NUTS), which is an adaptive form of a Hamiltonian Monte Carlo sampler, implemented through Stan 2.21.1 (Carpenter et al., 2017).
In Fig. 8, we compare our model results based on three of the previously introduced scenarios (I -III, see Table 1).The prior distribution of T plays an important role in the model, since longer durations will generally give higher likelihoods (as long as there is enough variation in the data to require it).Here, we have chosen to use a weakly informative prior distribution, T TMAX ∼ Beta(1, 2), that promotes a shorter duration without constraining the model more than necessary (see further discussion below).The results from Scenario I (Fig. 8A-B) show a maximum posterior marginal probability of T = 300 years.In addition, the results imply that durations of less than 100 years for the shield building phase are unlikely given the variations in the palaeomagnetic data.We note that this is true despite the model for Scenario I being unable to reproduce the shallow inclinations from flows 16-20.We note that the model does not penalise sequential outliers, which implies that this result is even less probable (or that the model is missing a component such as post-depositional tilting required to explain the data).In Scenario II (Fig. 8C-D), the more relaxed T MAX constraint means that the model is partially able to reproduce the shallow inclinations from flows 16-20.This leads to a bimodal posterior marginal distribution for T with high probabilities around both T ≈ 300 years and T ≈ 700 years.As with scenario I, the model for Scenario II shows almost zero probability for a total duration of less than 100 years.In Scenario III (Fig. 8E-F), the model is able to account for most of the variability in the palaeomagnetic data, as intended.Similar to scenario II, the relaxed constraints on T MAX leads to a bimodal posterior marginal distribution for T, with two peak probabilities around T ≈ 300 years and T ≈ 1250 years.
One of the main observations from the results presented in Fig. 8 is that the palaeomagnetic data appear to be incompatible with durations of T < 100 years, as suggested by Linnell et al. (2016).In the supplementary material we show that this observation is largely insensitive to the choice of prior distribution for T (Supplementary Material Fig. S3).Even when the prior distribution is defined so that there is a 95% prior probability that T < 100 years, the model still strongly suggests a longer duration for the shield building phase (90% probability T > 100 years), with a maximum posterior marginal probability of T ≈ 150 years.Without additional chronological constraints to those provided by the radiocarbon ages, we do not see any reasonable argument to propose stricter prior probabilities for T than this.
Another observation from Fig. 8 is that the shallow inclinations from flows 16-20 are incompatible with a short duration (T ≤ 505 years) implied by Scenario I.This implies that either: (i) the inclination data from flows 16-20 are systematically too shallow, perhaps due to an undetected reorientation of the lava flows after the NRM has been acquired, or that (ii) the youngest radiocarbon date found in the marine sediments below the analysed lava flows is too young, e.g., due to contamination with modern carbon.A third explanation is that the proposed upper bound of the shield-building phase (~550 years BP) may need to be revised, as suggested by Shane et al. (2013), i.e.-scenario IV.To test how compatible the results from the statistical model are with the various reference curves shown in Fig. 7, we construct preferred eruption chronologies for each scenario by calculating the mean posterior time passed between the individual flows Δt i .In Fig. 9, we compare the data and the model results plotted on the preferred chronologies for scenario I, III and IV to the different reference curves.As previously noted in section 5.2, the directional data appear to be more compatible with younger upper age of the shield building phase of 350 cal years BP, i.e.-scenario III and IV (Fig. 9B -C), as opposed to 550 cal years BP, scenario I (Fig. 9A).We also note that there is evidence for a brief period of shallower inclinations around 550 cal years BP, which is potentially compatible with both scenario I and IV.The observed variations in palaeointensity, although inherently more uncertain, appear to favour an upper age of 550 cal years BP, i.e.-scenarioI (Fig. 9D), but are also generally compatible with the longer eruption duration of scenario III (Fig. 9E).Overall, we conclude that, based on the differences between the reference curves and the associated model and data uncertainties, it is currently not possible to rule out any of the proposed scenarios based on the palaeomagnetic data.

Conclusions
The palaeomagnetic investigation of a sequence of lava flows from the western flank of Rangitoto has yielded 203 palaeoinclinations, 33 accepted absolute palaeointensity results and 39 accepted pseudo-Thellier palaeointensity results.Both the inclination and intensity data show significant, and often synchronous, variations between different lava flows, which suggest that these are related to changes in the past geomagnetic field.Attempts to constrain the duration of the volcanic activity in the shield-building phase based on comparisons with regional archaeo-/palaeomagnetic data and geomagnetic field model predictions are largely inconclusive due to large uncertainties in the reference data.We therefore constructed a statistical model to estimate the minimum duration required to account for the recorded field variations, based on prior information on the rate of change of the geomagnetic field.The results show that a short duration (<100 years) of the shield-building phase, as suggested by Linnell et al. (2016), is incompatible with the new palaeomagnetic data.The variations in the palaeomagnetic data alone would suggest a duration of 1000-1500 years (Fig. 8F), which is more consistent with predictions from the COV-LAKE model (Fig. 7).As indicated by Shane et al. (2013), some radiocarbon ages used to constrain the end of the eruption episode need to be re-evaluated, which could explain some of the discrepancy between the palaeomagnetic data and other age constraints.We also note that part of the variations in the palaeomagnetic data could have resulted from undetected re-orientation of the lava flows, e.g., due to rheological deformation, after the NRM was acquired.Based on the assumption that both palaeomagnetic and radiocarbon data are valid, our results suggest that the duration of the shield building phase was most likely in the range of 150-450 years (Fig. 8B).Such a long duration of activity at monogenetic volcanoes is unexpected and would be an additional factor for consideration in future hazard planning scenarios.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1. (A) Map of New Zealand.The red rectangle shows the area that is included in inset (C).(B): A map showing the location of the drill site (marked by a black diamond).Modified from Linnell et al. (2016).(C): A map of the Auckland region with relevant locations labelled.Areas covered with AVF deposits are also marked.

Fig. 3 .
Fig. 3. (A)Day plot for the samples analysed on the VFTB, colour-coded by compositional group(Linnell et al., 2016).Purple = geochemical group 1 (bottom of core), yellow = group 2, orange = group 3 and blue = group 4 (top of core).Day Plots(Day et al., 1977) plot the ratio of saturation remanent magnetisation to saturation magnetisation (Mrs/Ms) versus the ratio of coercivity of remanence to coercivity (Bcr/Bc).SD = single domain, PSD = pseudo-single domain and MD = multi domain.Theoretical SD/MD mixing curves (dashed lines) fromDunlop (2002aDunlop ( , 2002b) ) are shown for reference.Examples of hysteresis plots for (B) 99.4 m and (C) 122.0 m.The hysteresis curves were normalised to the maximum magnetisation (black lines) and corrected for paramagnetic influence (red lines).Temperature dependence of magnetisation for samples from (D) 57.2 m and (E) 122.0 m (experiments conducted in air), red = heating, blue = cooling.Examples of normalised susceptibility versus temperature for samples from (F) 121.05 m, (G) 3.73 m, (H) 41.5 m and (I) 34.0 m in an argon atmosphere, red = heating, blue = cooling.Temperature cycling experiments (heating curves as solid lines, cooling curves as dashed lines) for three identified magnetic phases: (J) between 90 • C and 170 • C (3.73 m), (K) between 320 • C and 460 • C (41.50 m) and (L) between 500 • C and 590 • C (3.73 m).
Fig. 3. (A)Day plot for the samples analysed on the VFTB, colour-coded by compositional group(Linnell et al., 2016).Purple = geochemical group 1 (bottom of core), yellow = group 2, orange = group 3 and blue = group 4 (top of core).Day Plots(Day et al., 1977) plot the ratio of saturation remanent magnetisation to saturation magnetisation (Mrs/Ms) versus the ratio of coercivity of remanence to coercivity (Bcr/Bc).SD = single domain, PSD = pseudo-single domain and MD = multi domain.Theoretical SD/MD mixing curves (dashed lines) fromDunlop (2002aDunlop ( , 2002b) ) are shown for reference.Examples of hysteresis plots for (B) 99.4 m and (C) 122.0 m.The hysteresis curves were normalised to the maximum magnetisation (black lines) and corrected for paramagnetic influence (red lines).Temperature dependence of magnetisation for samples from (D) 57.2 m and (E) 122.0 m (experiments conducted in air), red = heating, blue = cooling.Examples of normalised susceptibility versus temperature for samples from (F) 121.05 m, (G) 3.73 m, (H) 41.5 m and (I) 34.0 m in an argon atmosphere, red = heating, blue = cooling.Temperature cycling experiments (heating curves as solid lines, cooling curves as dashed lines) for three identified magnetic phases: (J) between 90 • C and 170 • C (3.73 m), (K) between 320 • C and 460 • C (41.50 m) and (L) between 500 • C and 590 • C (3.73 m).

Fig. 4 .
Fig. 4. Example of alternating field (AF) demagnetisation, thermal demagnetisation and microwave demagnetisation Zijderveld plots in core coordinates.Note that samples are oriented arbitrarily for thermal and microwave experiments and that for the AF demagnetisation experiments the declination is undefined.(A-D) and intensity decay plots (E-H).Examples of Arai plots (I-L) based on pseudo-Thellier, conventional Thellier and MWS methods.

Fig. 5 .
Fig.5.Downcore variations in rock magnetic and palaeomagnetic data alongside the lava flow stratigraphy (white = fill, grey = lava flow) and the geochemical zones (or compositional groups) as designated byLinnell et al. (2016).(A) Inclination (dashed line represents the GAD prediction), (B) Maximum Angular Deviation (MAD), (C) Angular deviation from the origin (AD), (D) The median destructive field of the NRM (MDF NRM , plotted in black) and the median acquisition field of ARM (MAF ARM , plotted in red), (E) Absolute Thellier (black) and calibrated pseudo-Thellier (red) palaeointensities (dashed line represents the intensity of the present geomagnetic field in Auckland, 53.7 μT) and (F) uncalibrated psuedo-Thellier palaeointensities (black dots), including rejected samples (black circles).

Fig. 6 .
Fig. 6.ΔR values surrounding the Hauraki Gulf embayment between 100 and 750 cal yr BP as compiled by Petchey and Schmid (2020).Positive ΔR values are displayed in bold and negative ΔR values are shown with standard font.The location of Rangitoto is marked with a red circle.

Fig. 7 .
Fig. 7. (A) Average inclination per lava flow (where 3 or more inclination values are available) with 95% uncertainty bounds based on the Student's t-distribution.N.B.not all flows are equal in size and some flows are more heavily represented in the material received for sampling.(B) Average absolute (black) and calibrated pseudo-Thellier (grey) palaeointensities per lava flow.(C) and D) independent inclination and intensity data from New Zealand, relocated to Rangitoto assuming a geocentric dipolar field: archaeomagnetic data (black) from Geomagia50.v3 with 2-sigma error bars, the Lake Mavora PSV record with 2-sigma uncertainty envelope (red) and NZPSV1k with 2-sigma uncertainty envelope (cyan).The archaeomagnetic data from Geomagia50.v3 without independent age control or with age uncertainties in excess of 2000 years were removed.The data fromTurner et al. (2020) were re-plotted using the independent ages provided in the original publication.The relative palaeointensity data from Lake Mavora were converted to absolute values using the scaling factor proposed byTurner et al. (2015b) prior to being relocated.Also shown are the predicted inclinations and intensities of two geomagnetic field models: pfm9k.2(blue) and COV-LAKE (green), with the mean model prediction shown as bold lines and 50 possible realizations as thin semi-transparent lines.The four scenarios (I -IV) for the duration of the shield building phase are represented by arrows above (C) and (D).

Fig. 8 .
Fig. 8.Results from the statistical model for scenarios I -III.The upper panels show the Rangitoto inclination data averaged by flow with 95% uncertainty error bars and the posterior distribution of the modelled inclination: mean model (black line) and the 95% credible interval (grey shaded area).The lower panels show the marginal posterior distribution of the total duration of the shield building phase T (grey bars) and the prior distribution (dashed black lines).In scenarios I and II, we used the probability density function of the calibrated ages from the respective radiocarbon dates using the Marine13 calibration curve with ΔR = − 91 ± 139 years (red line) are shown for reference.

Fig. 9 .
Fig. 9. Data and model outputs for scenarios I, III and IV plotted on their preferred eruption chronologies and compared to a selection of reference curves and global models.(A-C) The open points show the average inclination per lava flow (where 3 or more inclination values are available) and the error bars are 95% uncertainty bounds based on the Student's t-distribution.The black line is the mean result of the statistical model in each scenario, surrounded by the 95% credible interval (grey shaded area).(D-F) Absolute (black) and calibrated pseudo-Thellier (grey) palaeointensity results.The records/models shown are Lake Mavora PSV record with 2sigma uncertainty envelope (red), NZPSV1k with 2-sigma uncertainty envelope (cyan) and pfm9k.2(blue) and COV-LAKE (green), with the mean model prediction shown as bold lines and 50 possible realizations as thin semi-transparent lines.

Table 1
Summary of the dates that define the four scenarios discussed in this paper.