A new protocol using the viscous remanent magnetization (VRM) of boulders to date cataclysmic geological events such as tsunamis, glacial floods, and landslides is presented and its performance is assessed against two jökulhlaups (glacial floods) of known age in Iceland. High-intensity jökulhlaups have the ability to break off large boulders from bedrock and emplace and rotate them. These rocks originally carried a remanent magnetization parallel to the geomagnetic field during their formation. After being rotated by the flood, they acquire a VRM parallel with Earth’s magnetic field. In continuous thermal demagnetization experiments the unblocking temperature of the VRM can be determined, and subsequent rock magnetic VRM acquisition experiments can be used to establish a relationship between the unblocking temperature and the acquisition time, from which the time since the flood can be determined. The protocol was tested on 44 boulders from 2 historical jökulhlaups in Iceland and found to yield good order-of-magnitude estimates: 72 yr (confidence limits 11–360 yr) versus known 155 yr at the Sólheimajökull jökulhlaup and 290 yr (confidence limits 80–2300 yr) versus known 288 yr for the Kotarjökull jökulhlaup. The method can therefore be a valuable tool for future dating of cataclysmic events.

Floods with recurrence periods on historic time scales may pose an important natural hazard. These include not only storm floods and tsunamis, but also jökulhlaups, i.e., sudden high-volume glacial meltwater outbursts, which are common in Iceland, where they may have affected early settlement in medieval time (Smith and Dugmore, 2006). For risk management and planning, it is imperative to understand the history and recurrence of these and similar events. Dating methods exist, but have various shortcomings. Radiocarbon dating relies on the presence of organic material, cosmogenic radionuclide dating of flood deposits relies on fresh exposed surfaces (Icelandic rocks are mostly covered by snow), lichenometry relies on environmental conditions (e.g., air pollution), and tephrochronology relies on frequent close-by volcanic eruptions. Most rocks, however, contain small magnetic minerals that acquire a natural remanent magnetization (NRM) in the direction of Earth’s magnetic field at the time of their formation (Fig. 1). If it is subsequently reoriented during a flood, the NRM will in general no longer be aligned with the Earth’s magnetic field; the magnetic minerals then gradually acquire a new viscous remanent magnetization (VRM) in the direction of the ambient field that partially overprints the original NRM (Néel, 1949). For magnetically uniform, single-domain (SD) particles, the size of a VRM is a function of time, temperature, mineralogy, and grain-size distribution (Néel, 1949). Both the VRM and the original NRM can be recovered by demagnetizing the samples to progressively higher temperatures and measuring the remaining remanent magnetization vector: first, removing and identifying the VRM, and then the NRM. The temperature at which the VRM is completely removed is a function of VRM acquisition time and can be used to date the movement of the rock (e.g., Heller and Markert, 1973). Unlike other methods, VRM dating intrinsically reliant on the rocks, and is independent of external factors.

Until now, there was no generally accepted VRM dating method: Heller and Markert (1973) used alternating-field (AF) demagnetization to date the construction age of Hadrian’s Wall in northern England, Borradaile (1996) and Borradaile and Almqvist (2006) determined empirical calibration curves from dated archaeological material to date events of unknown age, and Kent (1985) and Smith and Verosub (1994) used a combination of stepwise thermal heating and room-temperature measurements to study limestone burial and landslides. Sato et al. (2014) applied the same technique to two tsunami-emplaced coral boulders: one yielded a paleomagnetic age of 1700 yr versus a radiocarbon age of 243 yr, and the other yielded 14 ka versus a maximum of 6 ka age of the formation of the reef. Muxworthy et al. (2015) measured the remanences of basalt boulders emplaced by floods in Iceland and Lake Bonneville (North America) at elevated temperatures as they were continuously heated and demagnetized. One estimate was accurate (80 yr versus a historically recorded 91 yr, Iceland), one was accurate to order of magnitude (15 ka versus a poorly constrained age estimate of 2.5 ka, Iceland), and one was incorrect (3.2 Ma versus 15 ka, Lake Bonneville), possibly due to large ambient temperature variations. Crider et al. (2015), using stepwise thermal demagnetization, could distinguish ages of four glacial moraines from Icicle Creek (Washington State, USA), with three ages in good agreement with cosmogenic ages (between 12 and 72 ka) and the oldest one (105 ka) being overestimated by VRM dating.

In this paper we develop a new robust protocol based on an objective and quantitative framework for VRM dating, and test its performance on two historical glacial floods (jökulhlaups) in Iceland.

A total of 44 boulders between ∼1 and 2 m diameter (Item DR1 in the GSA Data Repository1) were sampled from 2 locations in Iceland in August 2013 (Fig. 2), Sólheimajökull, and Kotarjökull. Sólheimajökull is a glacial tongue at the southern tip of Mýrdalsjökull that overlies Katla volcano, which erupted multiple times throughout recorded history, causing large-scale jökulhlaups, the last time in A.D. 1860 (Eliasson et al., 2006). At Kotarjökull, which is part of Vatnajökull, the last jökulhlaup occurred in 1727 (Thorarinsson, 1958). Jökulhlaups recur at irregular intervals, but we assume that the boulders that are still onshore were only moved during the last event (moreover, the most recent remagnetization likely completely or almost completely overprinted any previous VRM). The expected ages at the time of collection are 155 yr (Sólheimajökull) and 288 yr (Kotarjökull).

For each boulder, 5–8 independently oriented 1 cm cores were taken using a drill, and immediately stored in magnetically shielded containers until sample preparation and thermal demagnetization on Orion three-axis, high-temperature, low-field vibrating sample magnetometers (VSM) at Imperial College London and the Geomagnetic Observatory Borok, Russia. This instrument demagnetizes and measures strongly magnetic samples such as basalts of up to 1 cm in size, and was calibrated using a thermocouple cemented into a sample to ±1 °C accuracy. Hysteresis loops, backfield curves, and first-order reversal curves (FORC) were measured for sample characterization on high-field Princeton VSMs at Imperial College and at the Institute of Rock Magnetism, University of Minnesota, USA.

The principle of VRM dating is to relate the demagnetization temperature of the viscous remagnetization to the acquisition time post-flood. For SD particles, there is an expression to relate the temperature at the field location TA and the time tA of VRM acquisition to the demagnetization temperature TD in the laboratory experiment and the time scale tD of the experiment (Pullaiah et al. 1975). Assuming shape anisotropy dominating the magnetic remanence, this can be written
where TC is the Curie temperature of the magnetic mineral and τ0 is a material constant that is 10−8 to 10−12 s, but remains poorly constrained (Berndt et al., 2015); Sato et al. (2014) used 10−10 s and Muxworthy et al. (2015) used 10−9 s. Other works (e.g., Kent, 1985; Smith and Verosub, 1994) assessed a modified version of Equation 1 by Middleton and Schmidt (1982) that is now known to be inappropriate for VRM dating (see Item DR2). For VRM dating, one must determine all the parameters in Equation 1 and solve for tA, the age of the flood (or more generally, the redeposited material).

Demagnetization Temperature TD in Curved Demagnetization Plots

The demagnetization temperature TD is the unblocking temperature of the VRM, i.e., the inflection point in a demagnetization plot where the viscous remagnetization is fully removed. Often demagnetization plots do not show a clear single unblocking temperature, but can show significant curvature. In these cases, selecting the unblocking temperature visually as done by other VRM dating attempts is highly subjective (e.g., Muxworthy et al., 2015). Therefore, the selection of this point has been automated, following an approach similar to that in Crider et al. (2015): first, the demagnetization data are smoothed with a spline fit; second, the differential direction of the demagnetization vector is calculated; and third, the point of intermediate direction between the original magnetization and the viscous remagnetization is chosen as the unblocking temperature TD. This way, unblocking temperatures could be obtained even from strongly curved plots. We consider the point of intermediate direction the best choice for theoretical reasons outlined in Item DR3.

Demagnetization Time tD for Continuous Thermal Demagnetization

In conventional stepwise thermal demagnetization experiments, a sample is heated in zero-field to some temperature TD, and kept for a time tD (typically a few tens of minutes), after which the sample is cooled to room temperature again and its remaining remanence is measured. The grains with blocking temperatures below TD are thereby demagnetized and the process is repeated at successively higher temperatures. In this case the time tD can be directly inserted into Equation 1. This procedure, however, is not practical to use with temperature increments <10 °C, due to time intensiveness and instrumental accuracies. A difference in 10 °C in temperature, however, implies an order of magnitude in the age due to the logarithmic nature of the equation. Therefore, we used an Orion VSM capable of continuously heating in zero magnetic field (residual field < 100 nT), while continuously measuring the remanent magnetization, allowing for a 1 °C temperature resolution. As Equation 1 assumes a constant temperature over the time tD, a correction for the continuous heating is developed in Item DR4, yielding an effective time scale teff with
where r is the heating rate and W is the Lambert W function, which is defined as the solution of x = W(x)exp[W(x)].

Curie Temperature

Two rock magnetic quantities are required: (1) the Curie temperature TC, and (2) the attempt time τ0. The Curie temperature is easily determined by measuring thermomagnetic curves of the spontaneous magnetization MS(T) and determining the point of greatest curvature (Ade-Hall et al., 1965).

Viscosity Parameter: Effective Attempt Time t0,eff

Various approaches have been proposed to determine the attempt time τ0 (Berndt et al., 2015), yet it remains poorly constrained. As the age estimate is directly proportional to τ0, it is critical to determine it accurately. For this purpose, it is not the actual physical value of τ0 (the period between two successive thermal excitations) that is of interest, but rather an effective τ0,eff that accurately relates TD obtained from vector demagnetization plots to tA. These plots often include non-SD effects like multidomain (MD) behavior, magnetostatic interactions, and thermal alterations; therefore we have developed a method to empirically determine τ0,eff that when used with TD obtained from these plots best predicts the corresponding tA. Even though Equation 1 strictly only applies to ideal SD grains, using an effective τ0,eff partially corrects for errors introduced by pseudo-SD, MD, and interacting grains as the way τ0,eff is obtained realistically recreates the post-flood remagnetization process in the laboratory:

First, the sample is heated in the Orion VSM to above TC and cooled in a small-applied field similar to Earth’s magnetic field, creating a new thermoremanent magnetization (TRM). Second, the sample is thermally demagnetized by heating in zero field, continuously measuring the remanent magnetization MTRM(T). Third, the demagnetized sample is cooled in zero field to some temperature TA (between 100 and 300 °C), a small field applied, and the sample left in the field for a time tA (between 10 min and 1 day). During this procedure the sample acquires a VRM at known temperature TA and time tA. Fourth, the sample is again demagnetized, measuring the remanent magnetization MVRM(T) of the VRM. The last two steps were repeated (4–9 times) for different acquisition times tA and temperatures TA. For each of these experiments, a synthetic orthogonal-projection vector demagnetization plot was constructed plotting the TRM on one axis and the VRM on a perpendicular axis (Fig. 3). The demagnetization temperatures TD for the VRMs were determined using the algorithm described here. Equation 2 is used to obtain the demagnetization time tD from the known heating rate r. The quantities TA, tA, TD, and tD are then used in Equation 1 to solve for τ0 (as the experiment is repeated several times per sample, the τ0 value that minimizes the least-square errors is chosen). As the experiment approximates the natural TRM and post-flood VRM acquisition in the field, the natural and laboratory curvatures should be similar and the determined demagnetization temperatures comparable.

Field Temperature TA

The temperature TA in the field (post-flood) is taken from mean annual temperatures from climate data from A.D. 1961 to 2013 available from the Icelandic Meteorological Office (http://en.vedur.is/): 5.5 °C for the station Vík í Mýrda 20 km from Sólheimajökull, and 4.8 °C for the station Fagurhólsmýri 8 km from Kotarjökull.

Not all of samples are equally well suited for the dating. Rocks may have been realigned, but they may have also been affected by other events, e.g., weathering and lightning strikes. They may also contain a large variety of magnetic mineralogies with complex (un)blocking spectra from which it is difficult to recover the VRM. An objective set of criteria is needed to identify suitable samples and reject unsuitable ones, and to analyze the accuracy of resulting age estimates.

Directional Analysis

An emplaced boulder is expected to carry an original NRM in a random direction (high temperature) and a northward VRM (low temperature) (Fig. 1). Several samples were taken for each boulder, directions obtained from principal component analysis (Kirschvink, 1980), and the mean directions per boulder using Fisher (1953) statistics. The following criteria must be met.

  1. Northward trend of VRM: the direction of the VRM should tend northward; if not, then the VRM is likely a pre-flood remagnetization (Muxworthy et al., 2015). The direction rarely aligns perfectly with the north, because of (1) paleo–secular variation, (2) distortion of directions on continuous thermal demagnetization technique due to the temperature variation of MS(T), (3) the slow statistical process of VRM acquisition that tends to give less clear directions than TRMs, and (4) non-SD behavior (e.g., overlapping MD tails; Dunlop and Özdemir, 2000). Therefore, we accept boulders having mean VRM directions that are closer to the present-day geomagnetic north than their primary NRM direction.

  2. Clustering of VRM (Muxworthy et al., 2015): all the samples of one boulder should have roughly the same VRM direction, otherwise some samples may have altered, been subjected to elevated temperatures, or have complex mineralogies that do not reliably record a VRM (an α95 of 60° was used as a cutoff value).

  3. Clustering of original NRM: we introduce the new criterion that all independently oriented samples of one boulder should have roughly the same original NRM; otherwise, it would indicate that some samples were altered significantly (α95 of 60°).

Mineralogical Quality Criterion

Equation 1 is valid only for rocks containing a single type of magnetic mineral, but not for rocks containing a variety of different magnetic minerals, as those would acquire VRMs at different rates. A simple test is used to identify mineralogically suitable samples: the MS(T) plots should show a clear and unique Curie temperature. Samples with blurred out and/or exceptionally low TC (<300 °C) likely contain titanomagnetite assemblages of varying titanium content (our independent study), and were rejected.

Most samples showed two magnetic components (a VRM and an original NRM), either with a clear inflection point or with a strong curvature in the demagnetization experiments (Fig. 3, inset; Item DR5). In both cases, unblocking temperatures were obtained using the intermediate direction algorithm described here. Directional analysis found that most boulders carried a primary non-north clustering magnetization and a secondary northward clustering remagnetization, i.e., a VRM (Item DR5), but 9 of 44 boulders were rejected because their VRMs did not carry a secondary northward clustering remagnetization.

Most samples had a TC close to 580 °C, indicating magnetite, but 8 boulders were rejected for their blurred and low Curie temperatures, ∼200 °C (Item DR6). Hysteresis and FORC diagrams measured for most boulders generally indicated that samples from Sólheimajökull were more SD like, whereas those from Kotarjökull were more pseudo-SD and MD like, but no correlation between suitability for VRM dating and domain state was found (Item DR6). Viscosity experiments yielded median τ0,eff values of 5 × 10−8 s for Sólheimajökull (19 boulders) and 3 × 10−12 s for Kotarjökull (19 boulders), but showed no correlation with domain state (Fig. 3; Item DR7).

Using the τ0,eff value for each boulder in Equation 1 allowed for the calculation of post-flood acquisition times for each sample (Fig. 4). Of 20 boulders from Sólheimajökull, 12 passed the criteria outlined here; of 24 boulders from Kotarjökull, 14 passed (Item DR8).

Age Estimate and Statistical Error Analysis

Age estimates from individually dated samples yielded a large variance covering various orders of magnitude, including both overestimates and underestimates, which may be due to domain states, TA and TD uncertainties, or thermal or chemical alterations. While some of these ages are clearly incorrect (extreme values being 9 h and 1015 yr, i.e., larger than the age of the source rock or smaller than the time between sampling and measuring), no sample was rejected on the basis of the resulting age estimate, as this is the very quantity we aim to determine (doing so would introduce a sampling or confirmation bias); samples were selected purely on the basis of the selection criteria. How uncertainties propagate into the age estimate was investigated in Muxworthy et al. (2015); however, a full statistical treatment of sample variation is difficult because the variation is unlikely to be normally distributed: Equation 1 depends exponentially on TF and TL, but linearly on tL, and in a nontrivial way on the parameters used to calculate τ0. The underlying distribution is unknown; we therefore use a two-step bootstrap method similar to the one developed by Tarduno et al. (1990) to obtain both a flood-age estimate and uncertainty limits. The method (Item DR9) is based on (1) taking random resamples from the samples of each boulder to estimate the intersample variation for each boulder, and then (2) taking random resamples from all boulders to estimate the interboulder variation. From the resulting distribution, a median age and error limits corresponding to one standard deviation (i.e., the 16% and 84% quantiles) are obtained.

The final age estimate thus obtained is 72 yr for the Sólheimajökull flood, the actual known age of which is 155 yr, and 290 yr for the Kotarjökull flood, which has a known age of 288 yr. The confidence limits are 11–360 yr for the Sólheimajökull flood and 80–2300 yr for the Kotarjökull flood.

We have shown that our new protocol has the potential to successfully reconstruct the age of historic floods, with median ages closely approaching the real ages in the two test cases. It is more rigorous and has a more sophisticated error analysis than previous methods (e.g., Sato et al., 2014; Muxworthy et al., 2015). It is critical to have a sufficiently large sampling size because age estimates obtained from individual boulders can yield vastly different results; ∼20 boulders with 5 samples each is enough to obtain an order-of-magnitude estimate, but larger sampling sizes may reduce uncertainties. The method has been put on a sound theoretical foundation compared to previous studies, taking into account the heating-rate effect and correct choice of the rock magnetic parameters (effective attempt time, Curie temperature), controlling the VRM acquisition rate, and setting out a protocol for sampling, experimental procedure, data treatment, quality control, and error analysis. This should be a key step in establishing VRM dating as a tool for flood dating and subsequent risk assessment and mitigation.

We are grateful to V. Shcherbakov for use of the vibrating sample magnetometer at the Geomagnetic Observatory Borok, Russia, and for the visiting researcher fellowship to the Institute of Rock Magnetism (University of Minnesota), supported through the U.S. National Science Foundation and the Keck Foundation.

1GSA Data Repository item 2017098, additional theoretical derivations, detailed experimental data, and additional figures, is available online at http://www.geosociety.org/datarepository/2017/, or on request from [email protected].