Full‐Vector Paleosecular Variation Curve for the Azores: Enabling Reliable Paleomagnetic Dating for the Past 2 kyr

For archeomagnetic dating, high‐quality directional and intensity paleosecular variation curves are needed. The Azores Archipelago in the mid‐Atlantic Ocean provides a wealth of volcanic products erupted during the Holocene, making it an ideal location to (1) gather paleomagnetic data from well dated lava flows and (2) construct a paleosecular variation (PSV) curve that enables paleomagnetic dating of volcanic products with unknown age. Here, we present new full‐vector paleomagnetic data from Pico Island, and combine the new data with existing data from neighboring islands to construct a new full‐vector PSV curve for the Azores Archipelago. An extensive rock‐magnetic study underpins the quality of our paleomagnetic carriers. From Pico Island, we obtained 21 new mean site directions; and 15 paleointensity estimates with the multimethod paleointensity approach from 12 sites, the age was known for 14 and 10 sites, respectively. By bootstrapping the non‐Gaussian uncertainty estimates of the radiocarbon age calibrations and the confidence intervals associated with the direction and paleointensity estimates, we produce the first full‐vector PSV curve with confidence intervals for the Azores covering the past 2 kyr. The PSV curve reveals a period of low inclination between ∼900 and 1560 AD, with minimum values of 32°. The potential of our new full‐vector PSV curve is demonstrated by successfully dating five lava flows from Pico Island.


Geological Setting and Sampling
The nine volcanic islands of the Azores Archipelago (36. 9°N-39.8°N, 24.9°W-31.3°W) straddle the triple junction of the North-American, Eurasian, and Nubian lithospheric plates in the mid-North Atlantic Ocean (inset of Figure 1). Pico Island (447 km 2 ) has an elongated shape with a maximum length and width of 46 and 16 km, respectively. It is the youngest island of the Azores Archipelago (270 ± 150 ka; Demande et al., 1982) and is composed of three volcanic systems: Topo volcano, Planalto da Achada fissure zone and Pico volcano (Figure 1). Topo is an extinct volcano located in the southern part of the island that is deeply eroded and covered by younger lavas. Planalto da Achada fissure zone corresponds to a WNW-ESE-oriented ridge of scoria cones, ∼30-km long, that forms the eastern part of the island. Pico volcano is the youngest volcano of the archipelago (53 ± 5 ka; Costa et al., 2014) and constitutes the western part of the island. This central volcano is characterized by gentle lower slopes with numerous scoria cones, rising up to 2,351 m as a conical shape edifice with steep slopes toward the top (Nunes, 1999(Nunes, , 2020Zanon et al., 2020). The geological record of Pico Island reveals a high eruptive frequently during the Holocene (Nunes, 1999). Planalto da Achada fissure zone erupted at least 13 times in the last 2,000 years and Pico volcano erupted 22 times in the last 1,500 years (Nunes, 1999). After the island became inhabited in the 15th century, three historical eruptions have been recorded (Nunes, 1999;Zanon et al., 2020): Planalto da Achada fissure zone in 1562-1564 AD (Mistério da Praínha); Pico volcano in 1718 AD along two subaerial eruptive centers (Mistério de Santa Luzia on the north flank and Mistério de São João on the south flank) and a submarine vent offshore the south coast, and again in 1720 AD (Mistério da Silveira).
We sampled 22 sites from 21 independent cooling units on Pico Island ( Figure 1 and Table 1). For each site, ∼10 oriented samples were taken up to 20 m apart across each cooling unit to obtain a reliable paleodirection. For site PI07, the samples were taken closer together as the outcrop was rather small. Standard size paleomagnetic cores (2.5 cm in diameter and up to 15 cm in length) were collected using a petrol-powered BÉGUIN ET AL.  Table 2). The magnetization over temperature was measured on a horizontal translation Curie balance (Mullender et al., 1993). The Curie temperature was determined with the two tangent method and alternation of the sample was again determined by the nonreversibility of the magnetic signal.
Samples representative of seven sites show at least two dominant Curie temperatures, for sites PI05, PI10, and PI17 the Curie temperatures cannot be determined as the magnetization over temperature results in a gradual decay. The alteration temperature obtained by the magnetization or susceptibility vs. temperature experiments differ slightly for sites PI02, PI04, PI05, PI07, PI08, PI09, and PI21 (Table 2).
Based on the rock-magnetic behavior of the susceptibility and magnetization as function of temperature, the sites are divided in rock-magnetic groups after de   (Table 2). Rock-magnetic group L* is characterized by a sharp decrease in susceptibility from room temperature onward with a Hopkinson peak below room temperature. The remaining susceptibility at 150 °C is <80% of the room temperature susceptibility. Group H is the high-temperature group, it is characterized by an increasing or constant susceptibility from room temperature until at least 350 °C. Samples in group H retain 80% of their room temperature susceptibility at 400 °C, and have high Curie temperatures. Samples that show mixed behavior of both low and high Curie temperatures generally fall within group M, if the Hopkinson peak is below room temperature the samples are placed in group M*. Examples for each group are shown in Figure 2, results for all other sites are presented in supplementary information (S1 For all sites the follow information is provided: site number; location in UTM coordinates (UTM zone 26S); its age as provided (uncalibrated laboratory age for the 14C dating); the nature of the age constraint (14C: radiocarbon; historical: constrained based on historical observations); the median probability of the radiocarbon age (if applicable) after calibration with the INTCAL.13 curve; the one sigma probability intervals (68.3%) obtained by the radiocarbon calibration; and the reference dating site from Nunes (1999). Sites PI01 and PI02 are from the 1562-1564 historical lava flow, we use a mean age of 1563 AD for the eruption.  Dunlop (2002), with the exception of site PI04 with a high Bcr/Bc ratio and plots outside the PSD range and closer to the multidomain (MD) range.
Samples from sites with known age were checked for oxidation and exsolution. Carbon coated thin sections were examined with a JEOL JCM-6000 table-top scanning electron microscope (SEM) in backscatter mode, using an acceleration voltage of 15 kV. The scanned samples were categorized using the oxidation classes as described by Watkins and Haggerty (1967) (     Table 2 Rock-Magnetic, Paleodirection, and Paleointensity Results

Paleomagnetic Directions
Paleomagnetic directions were obtained from both thermal and alternating field demagnetization experiments. Four samples per site were thermally demagnetized in 13 temperature steps: 60, 100,150,200,250,300,350,400,450,500,530,560, and 600 °C, using an ASC TD48-SC furnace and were measured on a 2G cryogenic magnetometer (Figure 2, supporting information S4). For some samples, the magnetic moments at lower temperatures were beyond the dynamic range of the 2G cryogenic magnetometer; these low--temperature steps were not taken into account for the interpretation of the sample direction. Six to eight samples were demagnetized by increasingly higher alternating fields (AF) in 17 steps using a robotized 2G cryogenic magnetometer (Mullender et al., 2016); with peak fields of 2. 5, 5, 7.7, 10, 15, 20, 25, 30, 40, 50, 70, 80, 100, 150, 225, and 270 mT. All demagnetization data were analyzed using paleomagnetism.org (Koymans et al., 2016), site mean directions are calculated using Fisher statistics (Fisher, 1953), and the Vandamme cutoff (Vandamme, 1994) was used to identify potential outliers ( Table 2).
The obtained site means have small α 95 values (<7°), and the lowest obtained k value is for site PI02 with a value of 73.2. The directions of sites PI01 and PI02 should be the same since both sites are from the 1562 to 1564 AD historical lava flow: PI01 was sampled on the northern branch, whereas PI02 was sampled on the southern branch. The directions of these two sites are, however, not the same and this deviation needs to be explained. The two sites sampled on the 1718 AD historical lava flows, PI04 on the southern flank of Pico volcano and PI07 on the northern flank, have a common true mean direction (Koymans et al., 2020;Tauxe, 2010).

Paleointensities
Obtaining reliable and high-quality estimates of the past strength of the Earth's magnetic field is notoriously difficult. Recently, de Groot et al. (2013) showed the benefits of using multiple paleointensity methods in the same study. This multimethod approach combines classical (double) heating methods with isothermal paleointensity methods, minimizing the effects of thermal alteration, and dramatically increasing the success rate to obtain reliable paleointensity estimates from lavas. In this study, we applied three different paleointensity methods to our set of cooling units: thermal IZZI-Thellier (Tauxe & Staudigel, 2004;Yu et al., 2004), microwave Thellier (Hill & Shaw, 1999;Walton et al., 1993) reducing the amount of heat applied to the samples, and the calibrated pseudo-Thellier method (de Groot et al., 2013) performed at room temperature.
The principle for these three paleointensity methods is the same: the natural remanent magnetization (NRM) is compared to laboratory-induced magnetizations, these are imparted by known magnetic fields to assess the strength of the paleofield that imparted the NRM. The laboratory magnetizations are imparted thermally for the IZZI-Thellier technique, resulting in (partial) thermal remanent magnetizations ((p)TRMs). Microwave excitation is used to impart magnetization for the microwave Thellier technique, and for pseudo-Thellier strong asymmetric alternating magnetic fields impart anhysteretic remanent magnetizations (ARMs) in the sample.

Thermal Thellier
Based on the rock-magnetic analyses and the thermal demagnetization behavior, 12 out of 16 dated sites were chosen and initially five samples per site were subjected to the thermal IZZI-Thellier protocol ( Tauxe   The following rock-magnetic parameters are given per site: Curie temperature (Tc) and alteration temperature (Talt) obtained from magnetization-vs-temperature analyses, with Talt obtained from susceptibility-vs-temperature analyses provided within brackets, the room temperature (RT) susceptibility, the oxidation class (Ox. Cl.) according to ; Watkins & Haggerty (1967), and the rock magnetic group according to de Groot et al., (2015). The obtained paleomagnetic directions are given per site: the declination (dec), inclination (inc), 95 percent confidence interval (α95), precision parameter (k), and the number of samples measured (N), with the number of samples rejected by the Vandamme cut-off (rej.) between brackets. The obtained paleointensity results are given per protocol: Thermal Thellier, Microwave Thellier and pseudo-Thellier, the following results are presented per site: the number of samples that passed the selection criteria (n), number of samples measured (N), number of interpretations passing (different sets of) selection criteria (n_int), average paleointensity estimate from all passed interpretations (Paleoint.), standard deviation of the paleointensity estimate, and the standard deviation divided by paleointensity estimate. Paleointensity estimates that pass the site selection criteria are given in bold. Table 2 Continued including pTRM-checks (Coe, 1967), with a bias field of 60 µT. The samples were divided into two groups and the temperature steps were chosen accordingly: samples from rock-magnetic group L* were placed in the low-temperature group, with smaller temperature steps <400 °C and a total of 15 steps; the high-temperature group underwent additional steps >400 °C and went up to temperatures of 640 °C with a total of 16 steps. For three sites an additional 2-4 samples were measured in a second run, including pTRM-checks and a bias field of 40 µT. To ensure that the magnetic moments of the individual samples did not exceed the dynamic range of the 2G magnetometer, small samples (∼20 mg) were glued in borosilicate glass cups using quartz wool and KaSil glue. Samples were heated in an ASC TD48-SC thermal demagnetizer. The orientation of each sample was optimized so that the angle between the NRM and the laboratory field was always <90°, to suppress multidomain effects and help produce linear Arai plots (Paterson et al., 2015).
To discriminate between reliable and unreliable data, a multitude of selection criteria has been proposed.
We added the curvature criteria as defined by Paterson (2011) to all sets of selection criteria that did not already include this parameter, with |kʹ| < 0.35. To make the interpretation as objective as possible, the IZZI-Thellier data were analyzed with the auto-interpreter function of the Thellier-GUI (v. 3.13) (Shaar & Tauxe, 2013). All data were assessed using the six different sets of selection criteria (Figure 3), if a specific interpretation satisfied several selection criteria, it was taken into account multiple times, thereby this interpretation received more weight in the final result. The interpretations from samples of the same site were averaged, and the final site result was only accepted if the average paleointensity was calculated from at least three samples and the standard deviation divided by the mean paleointensity is below 20 % (bold values in Table 2).  (Table 2). For site PI06, two specimens from core PI06-X3 show very high paleointensity results, 94.8 ± 6.3 µT. The results from five samples from the same site but from different cores (PI06-X2 and PI06-X5) show paleointensity estimates of 44.1 ± 6.2 µT. Therefore, this site does not pass the site selection criteria.

Microwave Thellier
The measurement protocol in a microwave Thellier experiment is similar to the thermal Thellier technique: the original magnetization of the samples is stepwise replaced by a magnetization imparted in a known laboratory field. The technique used to magnetize and demagnetize the samples, however, is different. Microwaves are used to directly excite the magnetic spins in the samples, therefore the thermal energy to which the samples are exposed is relatively small compared to thermal Thellier techniques. This procedure may therefore reduce the effects of chemical alteration of the magnetic signals (Hill & Shaw, 1999). All measurements were conducted on the "Tristan" microwave system installed in the geomagnetic laboratory of the University of Liverpool (UK). In the microwave system, one single sample is processed before moving onto the next sample. This makes it possible to optimize the steps for each sample individually. First, one specimen from a site was stepwise demagnetized to obtain information on the demagnetization behavior, second, a sister specimen from the same site and core was subjected to an IZZI-Thellier paleointensity protocol, including pTRM-checks, with a bias field of 40 µT. To test the suitability of the samples to the microwave Thellier technique, samples from 10 different sites were subjected to a microwave Thellier experiment. For five of these sites the first, or first two, samples showed promising paleointensity estimates and additional samples were measured with a maximum of six samples per site. A total of 33 specimens were measured on the microwave system. The interpretation of the microwave data is done with the Thellier-GUI and the same six sets of selection criteria are used for the interpretation of the data. In total, 25 samples passed at least one of the selection criteria sets ( Table 2). The interpretations of some samples show ambiguous results, although passing the selection criteria. For example, sample PI07-06A ( Figure 3e) losses 80% of the NRM in the first four demagnetization steps, 68.1 W s, resulting in a paleointensity estimate of 185.9 µT. The quick demagnetization of this sample and the overprint on the NRM might indicate an effect of lightning although the other samples of this site do not show signs of lightning. As this is the only sample measured for this site, the behavior cannot be compared to other samples. Our microwave Thellier results highlight the importance of applying site selection criteria to the obtained individual sample interpretations. Of the 10 sites measured on the microwave system, two sites, PI01 and PI16, passed the site selection criteria.

pseudo Thellier
In a pseudo-Thellier experiment, samples are subjected to strong AFs, thereby avoiding thermally induced alteration caused by heating samples in conventional paleointensity methods altogether. While the thermal and microwave Thellier methods are absolute paleointensity methods, the pseudo-Thellier method is a relative paleointensity method (Tauxe et al., 1995). For lavas, absolute estimates of the paleofield can be obtained after calibration of the pseudo-Arai slope (de Groot et al., , 2013Lerner et al., 2017;Paterson et al., 2016). As the behavior of the pseudo-Thellier technique depends on the grain size distribution (Yu et al., 2003), samples suitable for calibration should satisfy the grain size selection proposed by de Groot (2013), which is the field at which half of the maximum ARM is gained for a sample, the B 1/2ARM value. In this study, we use the revised calibration formula from de Groot et al. (2016): B abs = 7.718 × |pTh-slope| + 14.600 µT, with 23 ≤ B 1/2ARM ≤ 63 mT. We use this calibration formula as it is based on natural remanent magnetizations rather that laboratory-induced thermal magnetizations as proposed by Paterson et al. (2016).
Samples from all sites were subjected to the pseudo-Thellier paleointensity experiment. Five to 11 specimens were first stepwise demagnetized using alternating fields up to 270 mT (Section 4). The same alternating field steps and a bias field of 40 µT are used to impart partial anhysteretic remanent magnetizations (pARMs). Lastly, the acquired ARM was demagnetized using the same field steps to check the linearity of the NRM-demag vs. ARM-demag. The data were analyzed with paleointensity.org (Béguin et al., 2020). The auto-interpretation function was used with a set of selection criteria specially designed for pseudo-Thellier experiments, pTh-SCRIT (Béguin et al., 2020). The focus of this set is on the linearity of the best-fit lines in all three plots (Arai, ARM-ARM, and Demag-Demag plots) and on the selection of the characteristic remnant magnetization (ChRM), where all signs of overprints on the NRM have to be removed to pass the selection. From the 146 measured specimens, 65 specimens passed the B 1/2ARM criteria. Of which seven sites hold realistic paleointensity estimates based on the pTh-SCRIT criteria set and the site selection criteria ( Figure 3 and Table 2).

Paleodirections
For all sites with oriented samples, a reliable paleodirection could be obtained ( Table 2). The samples were taken over a maximum distance of 20 m to average out local field anomalies and to obtain a more reliable paleodirection (e.g., Baag et al., 1995;de Groot & de Groot, 2019;Valet & Soler, 1999). This sampling strategy, however, might result in a lower k value than is expected for lavas. Nevertheless, our obtained k values are well above 70. Each cooling unit was sampled at one location with the exception of the 1562-1564 AD historical lava flow, which was sampled at the northern (PI01) and southern (PI02) branches. The paleodirections of PI01 and PI02, however, do not share a common true mean direction, meaning that these sites correspond to different cooling units or one of the sites might have moved postcooling, although there was no evidence for this in the field. Site PI13, with a mean calibrated radiocarbon age of 1543 AD, does show a common true mean direction with site PI02 (Koymans et al., 2020;Tauxe, 2010). Therefore, the direction of site PI02 is preferred over the direction of site PI01 for the 1562-1564 AD flow. The age of the PI01 flow is consequently considered unknown.

Multimethod Paleointensity Data From Pico Island
Seventeen sites were subjected to two or more paleointensity approaches, of which 12 sites passed at least one of the protocols and corresponding sample and site statistics; i.e., a site success rate of 70%. The three paleointensity methods rely on comparing the NRM of the specimen with a laboratory-induced magnetization, the technique used to induce the laboratory magnetization is different for each of the methods. The consistency of the paleointensity outcomes from the multimethod paleointensity technique can be considered as reliability check (e.g., Böhnel et al., 2009;de Groot et al., 2013;Monster et al., 2015). Site PI05 passed both pseudo-Thellier and thermal Thellier paleointensity protocols. Site PI16 passed all three paleointensity protocols. In both cases, the obtained paleointensity estimates from the different protocols are outside each other's confidence intervals. For these sites, the consistency between different methods is low. The thermal Thellier results of site PI16 is 24 µT higher than the obtained pseudo-Thellier. The microwave Thellier result for site PI16 (59.4 µT) is in-between the pseudo-Thellier and thermal Thellier result of 51.8 and 75.8 µT, respectively. For site PI05, the pseudo-Thellier result is 10 µT higher than the obtained thermal Thellier estimate. As all interpretations pass our selection criteria and are technically equally correct, it is impossible to imply which method yields the correct estimate of the paleofield. Sites with close ages, e.g., PI02 and PI13 or PI04 and PI07, show similar paleointensity estimates, within error, obtained by different methods; the thermal Thellier result of PI02 (44.0 ± 6.1 µT) and the pseudo-Thellier result of PI13 (41.9 ± 5.4 µT), also the pseudo-Thellier result of PI04 (39.3 ± 2.2 µT) and the thermal Thellier result of PI07 (43.2 ± 4.9 µT). This results in a high consistency between the paleointensity estimates obtained from different protocols for these sites.

Full-Vector Record for the Azores
The obtained paleomagnetic data from Pico Island can now be compared to the results of the nearby islands of São Miguel (Di Chiara et al., 2012, 2014b and Terceira , and existing geomagnetic field models ( Figure 4) (Jackson et al., 2000;Korte et al., 2011;Nilsson et al., 2014;Pavón-Carrasco et al., 2014;Thébault et al., 2015). For this comparison, all data are relocated (Noel & Batt, 1990) to a central point in the Azores Archipelago (38.5°N, 28°W), just northeast of the coast of Pico Island. Due to the maximum distance between the islands within the Archipelago the maximum relocation error is within 1° and 0.8 µT (Casas & Incoronato, 2007), although the expected relocation error is probably much lower. The directional data from the three islands show coherent results and coeval sites show similar directions within error. The overall trend of the data is covered by the geomagnetic models, although the high and low variations we find are obviously smoothed in the models. The data reveal a period of high declination between 600 and 800 AD and around 1600 AD which is not predicted by the geomagnetic field models. The most prominent mismatch of the geomagnetic field models and the data from the Azores is the inclination low (∼34°) between ∼900 and 1560 AD. This inclination low is obtained from seven sites from the islands of Pico and São Miguel spanning this time period (sites PI02, PI13, PI15, PI16, SM08, SM09, and SM10). The decline of the inclination could have started shortly after 780 AD, the first low inclination value sampled is from 934 AD (site PI16; 893-986 AD). This deviation from current geomagnetic field models amplifies the importance of the Azores Archipelago to constrain geomagnetic behavior in the mid-North Atlantic region.
The intensity record shows more scattered results, and the existing geomagnetic field models are generally within the data uncertainties with the exception of site PI19 around 1844 BC with the lowest paleointensity (31.2 ± 1.4 µT) from the Pico sites and the intensity high observed from two sites SM03 from São Miguel with a paleointensity of 93.0 ± 11.7 µT around 590 BC and PI14 with a paleointensity of 74.0 ± 4.6 µT around 870 BC. For these two high intensities, corresponding declination and inclination values agree well with the geomagnetic field models.

Paleosecular Variation at the Azores During the Past 2 kyr
Full-vector paleosecular variation (PSV) curves can be obtained at three levels: global, regional, and local. The global PSV curves can be obtained from archeomagnetic models e.g., CALS10k.   (Nilsson et al., 2014); CALS10k (Korte et al., 2011) including one standard deviation error envelops teal colored line; and the IGRF/gufm1 in black (Jackson et al., 2000;Thébault et al., 2015). The declination and inclination data from Pico are depicted as orange circles with the associated error bars calculated from the α 95 values (following calculations in Suttie and Nilsson (2019)). In the declination (a) and inclination (b) panels, data of Terceira are squares and from São Miguel are triangles; in the paleointensity panel (c), the shape of the data point reflects the paleointensity method used (see inset). The vertical error bars are the associated one standard deviation confidence intervals. The associated age errors are presented in Table 1 and supporting information S5.
used an existing global model, CALS7K.2 by Korte and Constable (2005), as an a priori model and used five published PSV curves to make a model for Europe to generate PSV curves at different locations. Regional PSV curves are generated by regional paleomagnetic models covering regions of continental scale, e.g., the regional model for the geomagnetic field in Europe of Pavón-Carrasco et al. (2008) using spherical cap harmonic analysis. To obtain local PSV curves, the paleomagnetic data for a small region or country can be interpolated to obtain the behavior of the paleofield for such small regions. For these local PSV curves, the paleomagnetic data must be relocated to a common point. Over the last decades different approached have been used to develop local PSV curves. Some curves are produced using hierarchical Bayesian modeling (e.g., Gómez-Paccard et al., 2012;Lanos, 2004;Lanos et al., 2005), while Thébault and Gallet (2010) introduced a bootstrap approach used by, e.g., Genevey et al. (2013) and Molina-Cardín et al. (2018). These approaches to derive local PSV curves have in common that they require a dense data input from which outliers can be identified and filtered from the data set.
The Azores Archipelago is located in a remote region in the middle of the North Atlantic Ocean, with the Canary Islands at ∼1,500 km and the Iberian Peninsula at ∼1,700 km distance. Since we find considerable deviations from existing geomagnetic field models, especially in our inclination data ( Figure 4) we restricted our data set to only include data from the Azores islands. All data with independent age constrains are included in our data set (supplementary information S5). Sites that are paleomagnetically dated in any of the previous Azores studies Di Chiara et al., 2012) were not considered. Site PI01 is not included in the data set as the paleomagnetic data of this flow shows inconsistent results when compared to other sites with coeval ages. As the data density is very low before 100 AD, we can only interpolate the past 2 kyr for our full-vector PSV curve. The available data are bootstrapped with a subsampling of 5,000 points within the uncertainty of the geomagnetic data and the (radiocarbon) age constrains. The calibration of radiocarbon age to calendar years is nonlinear as a result of changes in the rate of 14 C production in the upper atmosphere. The geomagnetic direction is parametrically sampled from the Fisher distribution of the site. The uncertainty of the intensity data is subsampled using a Gaussian probability function with two standard deviation. The age uncertainty is described by a probability density function as the distribution is non-Gaussian due to the radiocarbon calibration, obtained with the MatCal software (Lougheed & Obrochta, 2016). For each subsample set, we ran a cubic smoothing spline, and the final PSV curve results from the median of the 5,000 individual curves with 68% and 95% error envelopes ( Figure 5; supporting information S6). In addition to the Azores data set, synthetic data from the gufm1 model and the IGRF model were used to constrain the field expression from 1800 AD to the present-day field.
The resulting PSV curves ( Figure 5) describe the declination high around 600-700 and 1350-1650 AD and the inclination low around 900-1500 AD. For sites with coeval ages and differing paleomagnetic data, the smoothing spline fitted the average value, as all data was weighted the same. By repeating this average fit in the bootstrap, the estimated error envelop can result in a lower range than would be expected from the data. This is, e.g., the case for the inclination record around 1550 AD, where the error envelop is rather small compared to the data used around this age. When no paleomagnetic data are present the error envelops widen, e.g., in all three curves around 400-550 AD. The main observation from our new PSV curve is the sweep presented in the inclination record. This inclination low is not (profoundly) observed at the Canary Islands (e.g., de Groot et al., 2015;Kissel et al., 2015). Compared to European PSV curves (e.g., Batt et al., 2017;Molina-Cardín et al., 2018;Pavón-Carrasco et al., 2009) and the geomagnetic field models evaluated for the Azores Archipelago, again the inclination record stands out. The inclination record shows deviating behavior from 900 to ∼1700 AD. This shows the importance of a PSV curve in a remote location like the Azores Archipelago. The resulting PSV curves presented here can now be used to date volcanic products with unknown age that are within the time frame of the PSV curve, i.e., the past 2 kyr.

Paleomagnetic Dating with PSV Curves From the Azores
Seven sites (PI01, PI10, PI17, PI18, PI20, PI21, and PI22) were sampled from lavas with unknown ages. For these lava flows, an age range estimate is provided by the geological map of Nunes (1999) (Figure 1 and Table 3). The lava flows sampled at sites PI10, PI17, and PI18 are within the 2 kyr range of the PSV curve; while , and intensity (c) curves are shown with the 68% and 95% confidence error envelops. The data used to produce the curves are plotted in the appropriate graphs, symbols, and colors used are the same as for Figure 4; the thin colored lines present the geomagnetic field models for Pico Island. The non-Gaussian age error distributions are considered in the bootstrap. site PI21 might be within the range of the PSV curve but can also be up to 550 years older; and sites PI20 and PI22 fall outside the dating range of the PSV curve. We used the Matlab routine by Pavón-Carrasco et al. (2011) to estimate the age of the sites within the PSV dating range (supporting information S7), using our newly obtained PSV curves as master curve. The declination and inclination of the sites (Table 2) is compared to the PSV curves for the Azores. Sites PI01 and PI18 were the only flows that passed one of the paleointensity methods and for these flows the paleodirection together with the paleointensity (Table 2) was used to date the lava flows. The combined probability density functions of declination and inclination (and intensity for PI01 and PI18) result in multiple possible age ranges per flow (Table 3; supporting information S7).
Site PI01 was thought to be sampled on the northern branch of the 1562-1564 AD historical lava flow. However, the paleomagnetic data of this flow are inconsistent with site PI02 sampled on the southern branch of the same 1562-1564 AD flow and also inconsistent with site PI13, a radiocarbon dated flow with a calibrated age of 1543 AD. We thus assume that site PI01 was not sampled from the 1562 to 1564 AD lava flow but rather from an older flow below the historical flow. The paleomagnetic dating of this flow results in an age interval between 654 and 774 AD, validating this assumption. According to the geological map, the lava flows sampled at sites PI10, PI17, and PI18 are within the same age ranges (950-1450 AD). Site PI10 was sampled on the western part of the island in the lower slopes of Pico volcano, whereas sites PI17 and PI18 were sampled on the eastern part of the island in the Planalto da Achada fissure zone. The paleomagnetic dating of site PI10 result in a possible age between 1280 and 1440 AD with the highest probability at 1309 AD. The low inclination of PI10 fits the inclination low around this age interval. The obtained inclinations of sites PI17 and PI18 are too high to fit the inclination low of the Azores at this time interval, with measured inclinations of 59.8° and 57.0°, respectively. The ages of these flows are therefore not within the expected age range from the geological map. For site PI17, the paleomagnetic dating results in a single interval between 469 and 574 AD. The dating of site PI18 result in four possible intervals. Sites PI17 and PI18 should be within the same age range according to the geological map, from this we interpret that the age of PI18 is between 370 and 518 AD. Site PI21 was sampled on the western part of the island in the southern flank of Pico volcano, the age range provided by the geological map is −550-950 AD and this lava flow should be older than the lava flow sampled at site PI10. The measured inclination for site PI21 is 30.9°, which means that this flow must have been emplaced during the inclination low (∼900-1560 AD). The combined probability of the inclination and declination result in an age interval of 1184-1233 AD.

Conclusions
The volcanic products from Pico Island, in the Azores Archipelago, proved to be good recorders of geomagnetic field variations in the mid-North Atlantic region. We obtained paleointensity values for 12 of the 22 measured sites, of which 10 sites with known age. By combing our new data with existing data for the neighboring islands of São Miguel and Terceira we compiled the first full-vector PSV curve for the Azores region for the last 2,000 years. By using a bootstrap method, the non-Gaussian age distribution resulting from the calibration of radiocarbon ages was considered, together with the uncertainty of the paleomagnetic data. The data and PSV curves reveal an inclination low from ∼900 to 1560 AD, with minimum values of 32°. This inclination low is not previously constrained by global geomagnetic models and highlights the importance of a PSV curve for the Azores region. Furthermore, we illustrated the potential of paleomagnetic dating using our new full-vector PSV curve for the Azores by successfully constraining the ages of five lava flows from Pico Island. PI22 −3050-−550 --All flows are stratigraphically or historically constrained to the geological map age range. These broad ranges were narrowed to the paleomagnetic dating ranges provided; if the paleomagnetic dating yielded two age ranges, the one with the most probable age is given in bold. The most probable age is derived from the peak in the combined probability density function plot.