the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Direct U–Pb dating of carbonates from micronscale femtosecond laser ablation inductively coupled plasma mass spectrometry images using robust regression
Guilhem Hoareau
Fanny Claverie
Christophe Pecheyran
Christian Paroissin
PierreAlexandre Grignard
Geoffrey Motte
Olivier Chailan
JeanPierre Girard
Uranium–lead (U–Pb) dating of carbonates by laser ablation inductively coupled plasma mass spectrometry (LAICPMS) spot analysis is an increasingly used method in the field of geosciences, as it brings very strong constraints over the geological history of basins, faults or reservoirs. Most ages currently published are based on the measurement of U and Pb ratios on spot ablations, using nanosecond lasers coupled to sector field or multicollector ICPMS. Here, we test a new strategy for the U–Pb dating of carbonates from 2D isotopic ratio maps, based on the use of a robust regression approach in the data reduction workflow. The isotopic maps, with a minimum area of 0.65 mm^{2} (∼ 1000 pixels of 13×25 µm resolution), are obtained using a 257 nm femtosecond laser ablation system at a high repetition rate (500 Hz) coupled to a highresolution ICPMS. The maps commonly show significant variations in isotope ratios at the pixel scale, allowing the plotting of pixel U–Pb ratios in concordia or isochron diagrams and the calculation of U–Pb ages. Due to the absence of individual ratio uncertainties, the ages are calculated using MMrobust linear regression rather than the more commonly used Yorktype regression. The goodness of fit to the data is assessed by the calculation of the residual standard error (RSE) of the regression and by the calculation of a mean square of weight deviates (MSWD) on discretised data. Several examples are provided that compare the ages calculated by robust regression with those obtained by other techniques (e.g. isotope dilution, LAICPMS spot analyses and the pixelpooling approach). For most samples, characterised by high U concentrations (> 1 ppm), robust regression allows for the calculation of ages and uncertainties similar to those obtained with the other approaches. However, for samples with lower U concentrations (< 0.5 ppm), the ages obtained are up to 10 % too young due to pixels with high U ∕ Pb acting as leverage points for the regression. We conclude that the U–Pb ages calculated by the regression method tested here, although statistically robust, should be critically analysed before validation, especially for samples with low U concentrations.
 Article
(8806 KB)  Fulltext XML

Supplement
(25136 KB)  BibTeX
 EndNote
Since the seminal publication of Li et al. (2014), the uranium–lead (U–Pb) dating of carbonates by laser ablation inductively coupled plasma mass spectrometry (LAICPMS) has progressed significantly with applications as varied as the dating of tectonic fractures and veins (e.g. Roberts and Walker, 2016; Beaudoin et al., 2018; Parrish et al., 2018; Nuriel et al., 2019), carbonate deposition (Drost et al., 2018), speleothems (Woodhead et al., 2019) or cements filling the porosity of reservoirs (e.g. Li et al., 2014; Godeau et al., 2018). This is due, among other benefits, to the very good spatial resolution of the method, which allows one to detect and exploit submillimetrescale heterogeneities in U and the initial Pb concentration; to date several samples per day; and, due to the recent availability of carbonate standards, to correct matrixrelated offsets of U ∕ Pb ratios (Roberts et al., 2020). Given the ubiquitous nature of carbonates, the technique is likely to be associated with a growing number of studies in the future (Roberts et al., 2020).
Most carbonate dating LAICPMS studies are based on spot ablations of sizes ranging from 80 to 235 µm (e.g. Roberts and Walker, 2016; Parrish et al., 2018; Ring and Gerdes, 2016), where isochrons can be constructed from the combination of several tens of ablation craters made on the same crystal or, by default, on adjacent crystals when their size is too small (e.g. micrite). There is now enough hindsight on the applicability of the method to show that this approach allows precise ages to be obtained (< 3 % in the most favourable cases) (Roberts et al., 2020). In addition, the fact that most ablation cells can accommodate multiple samples simultaneously makes it possible to automatically analyse multiple samples per session, including primary and validation reference materials (RMs). However, as recently recalled by Roberts et al. (2020), the heterogeneous nature of carbonates requires prior petrographic and preferably geochemical characterisation of the samples to maximise the chances of successful dating. An initial screening of samples by LAICPMS to produce images of diagnostic trace element concentrations is particularly adapted to localise Urich and unaltered zones, in which precise and meaningful ages can be obtained (Roberts et al., 2020).
An alternative carbonate dating method based on age extraction from LAICPMS maps was recently presented by Drost et al. (2018). As before, the simultaneous analysis of many trace elements (including, of course, U and Pb isotopes) allows for the visualisation of the most favourable areas for dating. However, in this case, the concentrations of elements characteristic of detrital contamination (e.g. Zr and Rb) or those that allow one to define alteration zones or distinct calcite cement generations (e.g. U, Pb, Fe and Mn) are used to define pixel exclusion thresholds. These are used to filter out the values corresponding to microscopic inclusions of other minerals (e.g. clays) and to altered zones, and/or to isolate the cement generation targeted for dating. The isotopic U and Pb values of the remaining pixels are used to build isotopic ratio maps, from which a U–Pb age can be calculated. The methodology of age determination proposed by these authors is based on the sorting of pixel ratio values used in dating (^{238}U ∕ ^{206}Pb–^{207}Pb ∕ ^{206}Pb for a TeraWasserburg, TW, diagram) based on ^{238}U ∕ ^{208}Pb and ^{235}U ∕ ^{207}Pb pixel ratios and their pooling (i.e. discretisation) to subsets (“pseudoanalyses”) of a similar number of pixels for which the mean values and associated standard errors (i.e. standard deviations of the mean) are calculated. These values are then plotted in TW, 86TW (^{208}Pb_{common} ∕ ^{206}Pb versus ^{238}U ∕ ^{206}Pb; Parrish et al., 2018) and isochron diagrams (^{206}Pb ∕ ^{208}Pb_{common} versus ^{238}U ∕ ^{208}Pb_{common}) to calculate the age using Yorktype regressions as implemented in Isoplot (Ludwig, 2012) or IsoplotR (Vermeesch, 2018).
In this contribution, we test the pros and cons of an alternative approach for the calculation of U–Pb carbonate ages from isotopic maps. Rather than the pooling of pixel values proposed by Drost et al. (2018), the pixel isotopic ratios of interest are directly plotted in the concordia or isochron diagrams, and a statistically robust linear regression that is extremely robust to outliers is processed through the pixels to obtain the age. The visualisation of all pixel ratio values in the TW or isochron diagrams allows the user to visually assess the quality of data, guided by a more common statistical parameter such as the residual standard error (RSE). Discretisation is also used to calculate pseudoMSWD (MSWD refers to the mean square of weight deviates) values that permit the assessment of the linearity of pixel data. The robust regression method appears to be ideally suited to isotopic maps obtained after pixel filtering based on trace element concentrations, as proposed by Drost et al. (2018), or if this is not feasible, using pixel colocalisation, as proposed here.
2.1 Choice of samples
Six samples of already dated lacustrine and reefal carbonates, and of tectonic veins, have been chosen to test the usefulness of the robust regression approach. Among these, four samples were analysed by the LAICPMS isotope mapping approach as part of this study. These samples are as follows:
 i.
a lacustrine limestone (Long Point; Duff Brown Tank locality in the Colorado Plateau, USA), precisely dated by Hill et al. (2016) at 64.0 ± 0.7 Ma (2 s) by U–Pb methods using isotope dilution (ID) multicollector inductively coupled plasma mass spectrometry (MCICPMS), and labelled Duff Brown in the following. This sample serves as a validation RM in the analytical procedure and has also been demonstrated to be suitable for U–Pb dating using image maps by Drost et al. (2018);
 ii.
a tectonic calcite vein (BH14) from the Bighorn Basin (Wyoming, USA) dated by U–Pb LAICPMS spot analyses (i.e. range of spot ablations) at 63.0 ± 2.2 Ma (MSWD = 1.6) by Beaudoin et al. (2018) and at 61.2 ± 2.9 Ma (MSWD = 4.1) in our laboratory (the detailed methodology is presented in Fig. S1 in the Supplement), with the same calcite RMs as the present study (WC1 and Duff Brown);
 iii.
a reefal limestone sample from the southern Pyrenees (PXG201), located in the upper part of the Pamplona marls in the Atarés anticline. This sample has been dated at 40.5 ± 2.3 Ma (MSWD = 0.9) in our laboratory by LAICPMS spot analysis. A micropaleontologic study from Canudo and Molina (1988) performed close to the sampling area suggests an upper Bartonian to Priabonian age for the deposits (∼ 38.5–36.0 Ma), whereas the paleomagnetic study of Oms et al. (2003) is more consistent with a Bartonian age (∼ 40.0–38.5 Ma);
 iv.
a lacustrine carbonate from the southern Pyrenees (PXG322), which is part of the first continental deposits of the Campodarbe Formation flanking the western limb of the Boltaña anticline (eastern Jaca Basin). Previous magnetostratigraphic and tectonosedimentary studies point to a Bartonian–Priabonian age of deposition (Mochales et al., 2016). A similar age was obtained by U–Pb LAICPMS spot analysis in our laboratory, although with poor statistics (37.0 ± 3.0 Ma; MSWD = 6.3).
The two remaining samples were previously analysed by Drost et al. (2018) and Roberts et al. (2020). The raw pixel isotopic ratios retrieved from their LAICPMS isotopic maps (Kerstin Drost, personal communication, 2020) were processed with the robust regression approach. These are
 v.
a calcite spar from a Carboniferous–Permian lacustrine limestone (sample JS4, run JS41), dated by Drost et al. (2018) at 300.8 ± 3.7, 299.6 ± 4.0 and 298.9 ± 5.6 Ma with their pixelpooling approach (TW, 86TW and isochron diagrams respectively). This age agrees with several depositional age constraints (see Drost et al., 2018 for details).
 vi.
a tectonic vein from the Bighorn Basin (Wyoming, USA) dated by U–Pb LAICPMS spot analyses at 59.5 ± 2.7 Ma (MSWD = 2.5) by Beaudoin et al. (2018). Roberts et al. (2020) obtained a similar age of 61.0 ± 1.7 Ma with the pixelpooling approach of Drost et al. (2018).
Finally, three additional, undated samples of calcite cements found in tectonic breccias and veins affecting Tithonian limestones of the northern Pyrenees (France) are also presented (ETC2, ARB and C6265D5). They provide examples of the application of the technique to the resolution of the diagenetic evolution of areas with complex tectonic evolution.
2.2 Analytical strategy
2.2.1 LAICPMS instrumentation
All the samples were analysed with a femtosecond laser ablation system (Lambda3, NEXEYA, Bordeaux, France) coupled to a sector field inductively coupled plasma mass spectrometer (SFICPMS) Element XR (Thermo Fisher Scientific, Bremen, Germany) fitted with the Jet Interface. The laser is fitted with a diodepumped Yb:KGW (potassium gadolinium tungstate) crystal laser source (HP2, Amplitude Systèmes, Pessac, France). The pulse duration is less than 400 fs at 257 nm. The laser source can operate within a wide range of repetition rates (1 Hz to 100 kHz) and energy ranging from 200 µJ per pulse below 1 kHz to 1 µJ at 100 kHz at 257 nm. Complex trajectories can be realised by moving the laser beam (15 µm diameter at full energy) across the surface of the sample using the rapid movement of galvanometric scanners combined with a high repetition rate (Fig. 1a; Ballihaut et al., 2007, Aramendia et al., 2015).
The aerosol produced by the ablation was carried to the ICPMS by a tube ($\mathrm{1}/{\mathrm{16}}^{\prime \prime}$ internal diameter) using a helium stream (600 mL min^{−1}), with the exception of the first session in October 2018 where argon (Ar) was used (650–700 mL min^{−1}). Measured wash out time of the ablation cell was ∼ 500 ms for helium gas and ∼ 1 s for Ar gas considering the 99 % criterion. To improve sensitivity, 10 mL min^{−1} of nitrogen was added to the twister spray chamber of the ICPMS via a tangential inlet while helium flow was introduced via another tangential inlet located at the very top of the spray chamber. Measurements were performed under dry plasma conditions. The femtosecond laser ablation inductively coupled plasma mass spectrometry (fsLAICPMS) coupling was tuned on a daily basis in order to achieve the best compromise in terms of sensitivity, accuracy, particle atomisation efficiency and stability. The additional Ar carrier gas flow rate, torch position and power were adjusted so that the U∕Th ratio was close to 1±0.05 when ablating the NIST SRM 612 glass. Detector crosscalibration and mass bias calibration were checked daily using the appropriate sequence of the Element software. The laser and ICPMS parameters used for U–Pb dating are detailed in Appendix A. For sample screening (first type of isotopic map, see Sect. 2.2.2), only ^{238}U, ^{207}Pb and ^{206}Pb were selected, reaching a total mass sweep time of about ∼ 500 ms. For the maps used for dating, ^{238}U, ^{232}Th, ^{208}Pb, ^{207}Pb and ^{206}Pb were selected, reaching a total mass sweep time of about ∼ 60 ms. The analysis of trace elements with lower atomic masses, for the purpose of pixel selection as done by Drost et al. (2018), was not attempted due to the inability of the SFICPMS to achieve reasonable total mass sweep times (due to magnet mass change). Based on NIST612, the limit of detection as of October 2019 was ∼ 1.3 ppb and ∼ 0.2 ppb for ^{206}Pb and ^{238}U respectively. The limit of quantification was ∼ 4.3 and ∼ 0.8 ppb for ^{206}Pb and ^{238}U respectively.
2.2.2 Analytical workflow
Two types of isotopic maps were produced:
 i.
The first one (screening map) was designed to identify the most favourable areas for dating (high U ∕ Pb), over a large surface of the sample. In this configuration, samples were first ablated along lines of 6.5 to 8.1 mm width at a repetition rate of 300 Hz, without precleaning. The lines of 50 µm height were obtained using a back and forth movement of the laser (at 10 mm s^{−1}) combined with a stage movement rate of 100 µm s^{−1} (Fig. 1a), corresponding to 64.5 to 81 s of analysis per line, followed by a 15 s break between the lines (to allow for ICPMS data processing). The lines were separated by a distance of 100 µm, leaving a gap of 50 µm between them. The number of lines was 40 to 56, resulting in a duration of 53 to 78 min for a complete map of a surface of between 26 and 38.1 mm^{2}. The maps produced (^{238}U ∕ ^{206}Pb) are presented in Fig. S2.
 ii.
The second one (dating map) was designed to date the samples. Therefore, they were performed in a selected area (from screening map and after being repolished) and consisted of linear scans of 0.8 to 1.6 mm width at a repetition rate of 500 Hz. These lines of 25 µm height, separated by a distance of 25 µm so that they are immediately adjacent to each other, were obtained using a back and forth movement of the laser (at 5 mm s^{−1}) combined with a stage movement rate of 25 µm s^{−1}, corresponding to 32 to 64 s of analysis per linear scan, followed by a 15 s break. The number of lines was 25 to 27, resulting in a total analysis time ranging from 19 to 38 min for a complete map of a surface of between 0.5 and 1.1 mm^{2} (Fig. 1b). Before analysis, the samples were precleaned with the laser using a stage movement rate of 200 µm s^{−1}. The unknowns were bracketed with the NIST SRM 612 for lead ratios and by the commonly used WC1 calcite RM (age 254.4 ± 6.4 Ma; Roberts et al., 2017) for the Pb ∕ U ratio, using the standardisation method of Roberts et al. (2017). Two unknowns could be analysed in a row between the RMs. The primary and validation RMs were analysed under conditions similar to the unknowns, except that the isotopic maps were of a smaller surface (∼ 0.2 mm^{2}), corresponding to an analysis time of ∼ 5 min. The ablation depth was between ∼ 25 and ∼ 35 µm.
2.3 Data processing and uncertainty propagation workflow
2.3.1 Construction of isotopic maps and blank correction
For each map, the first line corresponds to the ICPMS signal recorded without laser ablation in order to calculate average blank values for the different masses. The matrices of raw isotopic values (in counts per seconds) were then built using a homemade Python code including blank correction (subtraction of average blank for each data point). Visualisation as maps were processed with the Fiji freeware. All isotopic maps are imported as an image stack and appropriately cropped. To consider the wash out time, the number of pixels was averaged by eight along each linear scan in order to reach a value of one pixel per 500 and 1000 ms for He and Ar respectively. The isotopic ratio matrices were then calculated.
2.3.2 Normalisation of the Pb ∕ Pb ratios and drift correction of the U ∕ Pb isotopic ratios
For the NIST SRM 612 glass, Pb and U isotopic ratios were calculated as the robust mean of the ratios over the entire image, using the “Huber” function of the “Statsmodels” Python library. The robust mean is obtained through an iterative process that uses a Huber loss function to strongly reduce the weight of any outlier (i.e. spike) onto the final result, without any manual outlier rejection. The standard error (i.e. standard deviation of the mean) for each ratio of the NIST SRM 612 glass was also calculated from the standard deviation which is an output of the function. For pixel Pb ∕ Pb ratios of the unknown and WC1, normalisation (correction of mass bias) was based on a standard bracketing approach, using the NIST SRM 612 analysed immediately before and after the unknown. The reference values were taken from Woodhead and Hergt (2001). For U ∕ Pb ratios, correction for drift was based on a similar approach (standard bracketing using NIST SRM 612 glass). For WC1 calcite, an age was calculated from the combination of both isotopic maps that bracket the unknown sample. The age, obtained with the robust regression method described in detail below, was used to calculate the correction factors for the U ∕ Pb ratios, following the approach of Roberts et al. (2017). The correction factor was calculated for a TW regression, with the ^{207}Pb ∕ ^{206}Pb anchored to a value of 0.85 as determined from isotope dilution methods (Roberts et al., 2017). The U ∕ Pb isotopic ratios of all pixels of the unknown sample isotopic maps were corrected by the obtained correction factor.
2.3.3 Age calculation of the unknown
All regressions were performed using a statistically robust linear regression approach based on a Tukey biweight loss function. All of the pixel isotopic values were first plotted in the appropriate plot (TW, 86TW and isochron plots). For the 86TW and isochron plots, the amount of common ^{208}Pb (^{208}Pb_{c}) must be known to calculate the ^{208}Pb_{c} ∕ ^{206}Pb, ^{206}Pb ∕ ^{208}Pb_{c} and ^{238}U ∕ ^{208}Pb_{c} ratios. For each pixel, the ^{208}Pb_{c} was calculated from the amount of radiogenic ^{208}Pb obtained from measured ^{232}Th and the age of precipitation given by the TW plot. This heuristic approach is similar to that proposed by Parrish et al. (2018). The regression was then made using the “lmrob” package in R (Maechler et al., 2019). In the chosen setting (KS2014) the package allows one to automatically perform an MMrobust regression, which has a high breakdown efficiency and is demonstrated to be highly efficient against outliers (Koller and Stahel, 2011, 2017). Based on an iterative process, the method automatically attributes a weight to each point, with the lowest weights corresponding to the outliers. The final regression is then a weighted linear fit through the points. Finally, the age calculation used the classical geochronological equations, corresponding to the lower intercept of the regression line with the concordia in the TW space (with the x axis in the TW86 space) and to the slope of the regression line in the isochron diagram.
2.3.4 Uncertainties
In the commonly used Yorktype weighted linear regressions, the weight of each point is inversely proportional to its uncertainty (i.e. size of the ellipse), where the latter commonly corresponds to the quadratic addition of analytical uncertainties (Horstwood et al., 2016). In contrast, the pixel isotopic ratios used for the robust regression have no uncertainty, as they are directly calculated from the number of counts of the relevant masses. Therefore, all uncertainties (analytical or systematic) were added quadratically to the age uncertainty obtained from the confidence interval of the robust regression. Following the recommendations of Horstwood et al. (2016), these uncertainties were added to the U ∕ Pb ratio uncertainty corresponding to the age, rather than to the age uncertainty itself. The analytical uncertainties considered are those calculated for the homogeneous primary reference material (NIST SRM 612 glass). Due to the standard bracketing approach used (see above), no excess variance could be calculated. The systematic uncertainties are the decay constant uncertainty of ^{238}U (0.1 %, 2 s) and the ^{238}U ∕ ^{206}Pb ratio uncertainty of WC1, as estimated by Roberts et al. (2017) (2.7 %, 2 s), and a longterm excess variance taken as 2.0 % (2 s). Due to the low number of sessions carried out with this approach, the longterm variance could not be reliably calculated. This value was therefore chosen by analogy with other studies (e.g. Horstwood et al., 2016; Beaudoin et al., 2018; Guillong et al., 2020), but it is likely to change in the future.
2.3.5 Estimation of the goodness of fit
In geochronology, the parameter usually most used in the estimation of goodness of fit is the mean square of weight deviates (MSWD), or reduced chisquare, whose value close to one is interpreted as a strong indication of the quality of the regression. The calculation of the MSWD can, by definition, only be performed from values associated with uncertainties, which is not the case for the points used for the robust regression. Two steps have been adopted in this work to quickly assess the goodness of fit of the robust regressions from statistical parameters. The first one is the calculation of the F test of overall significance. If the F test is higher than the significance level (chosen at 0.05), the regression is rejected and the dating is considered unsuccessful. If this first step is successful, the second step considers several parameters that help with assessing the quality of the regression. One parameter is the residual standard error (RSE) of the robust regression, an output of the lmrob method, which provides an estimate of the mean dispersion of the points around the regression in absolute value of the y axis. The smaller the RSE, the smaller the dispersion. Based on our experiments with analysed carbonate samples, successful samples have an RSE lower than ∼ 12 % of the y intercept (Pb_{0}) value (i.e. < 0.10 for a common Pb value of 0.85 in TW plots). The second step is the calculation of a MSWD from a discretisation of the points used for regression, which we call MSWDd. For each type of plot, the points are discretised into sets along the regression line obtained with the robust regression. For each set, the number of points always corresponds to the equivalent of 30 s of analysis (i.e. 60 pixels for He and 30 pixels for Ar), which is the most common duration of laser ablations in spot LAICPMS studies. The weighted mean and standard error of the sets are calculated (i.e. ellipses). A MSWDd value close to one is expected for sets well aligned along the regression. Based on our experiments, the ages obtained with MSWDd values above ∼ 2.5 should be questioned. This approach is close to the one used by Drost et al. (2018) for their age calculations, except that (i) discretisation is not performed using the ^{235}U ∕ ^{207}Pb or ^{238}U ∕ ^{208}Pb ratios, (ii) the MSWD is not based on the use of Yorktype regression and (iii) the number of pixels per ellipse is not adjusted opportunely for each sample to tend toward a better MSWD value. In parallel with the calculation of the MSWDd, a TW plot with the corresponding ellipses, as well as a representation of the running mean with 60 and 30 pixels for He and Ar respectively, is also systematically produced to assist in the user's assessment of the quality of the point alignment and, thus, the validity of the linear regression.
Two additional checks are performed to ensure the validity of the results. First, the ages obtained according to the three types of diagrams (TW, 86TW and isochron) must be identical within uncertainty. Second, an age is calculated from the subdivision of the isotope maps into squares that we call “pseudospots” and that mimic conventional LAICPMS spot analyses. In detail, they have dimensions of 7×6 pixels, equivalent to 175 µm × 150 µm and 21 s of analysis, except for sample BH14 (October 2018) where they have a dimension of 4×6 pixels (100 µm × 150 µm; corresponding to 24 s). The mean and standard error of the ^{238}U ∕ ^{206}Pb and ^{207}Pb ∕ ^{206}Pb ratios are calculated from the pixels of each pseudospot, and the age is obtained in the TW space using IsoplotR. The age calculated with the robust regression must be, within the uncertainties, identical to that obtained from the pseudospots.
2.3.6 Colocalisation study
The lack of trace element concentration data prevents the adoption of a pixel selection approach as presented by Drost et al. (2018). In order to select the most favourable areas of the isotope maps to obtain precise ages, a pixel value colocalisation approach using the ImageJ plugin ScatterJn (Zeitvogel and Obst, 2016) was used. The pixel values of the ^{207}Pb ∕ ^{206}Pb and ^{238}U ∕ ^{206}Pb isotopic maps were plotted in a TW scatter plot, and areas selected in the spatial domain (maps) were highlighted in the plot (and vice versa). For all samples analysed, the spatial location of (i) outliers highlighted by the robust regression approach and (ii) visually aligned high ^{238}U ∕ ^{206}Pb ratio points in the TW plot were identified. In either case (or both), if the identified pixels are located in close spatial proximity on the map, a subset of the map can be selected and the age can be calculated from the pixels in that subset. It should be noted, however, that due to the small number of pixels selected, the age obtained is not necessarily more precise than using the entire map (as also highlighted by Drost et al., 2018, regarding their approach).
2.3.7 Comparison with the method of Drost et al. (2018)
In order to compare the efficiency of the robust regression method with that of Drost et al. (2018), we also carried out a processing of isotope ratio data directly inspired by these authors. It consists of sorting the pixel ratio values using the ^{207}Pb ∕ ^{235}U ratio for preCenozoic samples or the ^{238}U ∕ ^{208}Pb ratio for Cenozoic samples, clustering the data in discrete steps of a given number of pixels (here 60 pixels, or 30 pixels for sample BH14, corresponding to 30 s of signal), before calculating the mean and its uncertainty for each cluster. The age is then calculated in TW, 86TW and isochron diagrams using IsoplotR. The stated age uncertainties include the same systematic uncertainties as for the robust regression.
For all samples considered in this study, the three types of diagrams used for age calculations by the robust regression (TW, 86TW and isochron), the TW plots resulting from the use of the pixelpooling approach Drost et al. (2018) method and the plots produced from image discretisation (pseudospots), are presented in the Supplement (Fig. S3). The ages calculated by the different approaches are also listed in Table 1. Finally, the pixel values of the isotopic images for the masses 238, 232, 208, 207 and 206 are presented in Table S1.
3.1 Age determination on previously dated samples
3.1.1 Duff Brown
The ages were obtained on two images chosen as an example. The ages calculated from pseudospots in the TW space are imprecise, with values of 70.7 ± 20.2 Ma and 54.8 ± 13.3 Ma (MSWD of 4.4 and 2.8 respectively; Fig. S3; Table 1). This is due to the lack of spread of the data points (i.e. ellipses) along a linear array, which results in poorly defined slope and intercepts of the regression line. Indeed, anchoring the Pb_{0} to the value calculated from the data of Hill et al. (2016) (0.738 ± 0.020, 2 s) leads to more precise ages of 63.6 ± 2.5 Ma (MSWD = 4.2) and 66.3 ± 2.9 Ma (MSWD = 4.0), similar to the reference value within uncertainties (Fig. S3). The same observations can be proposed for the robust regression method. The ages found without anchoring the common lead value are 54.9 ± 2.5 Ma (Fig. 2b, d) and 57.4 ± 2.4 Ma (Fig. S3) (MSWDd of 8.5 and 4.1 respectively), whereas using the fixed Pb_{0} value leads to ages of 63.1 ± 2.1 Ma (Fig. 2c) and 65.9 ± 2.2 Ma (Fig. S3), similar to the reference value within uncertainties. Note that the maps used for age calculation (Fig. 2a) are smaller (0.21 mm^{2}, ∼ 630 pixels) than those made for other samples, as Duff Brown was used as a validation RM during the tests. In this case, the very high MSWDd values (39.1 and 29.5) are due to the low ^{207}Pb ∕ ^{206}Pb values of the pixels the closest to the y axis (i.e. with lower ^{238}U ∕ ^{206}Pb value) compared with those expected from the anchored regression. Finally, the ages obtained with the method of Drost et al. (2018) are identical to those from the robust regression, within the limits of uncertainties. The ages obtained without anchoring the Pb_{0} value are 61.9 ± 3.4 Ma (MSWD = 3.1; Fig. 2e) and 61.2 ± 2.5 Ma (MSWD = 1.3; Fig. S3), which are close to the expected value. By anchoring the Pb_{0} value, the ages are of 63.9 ± 2.2 Ma (MSWD = 4.5) and 65.8 ± 1.0 Ma (MSWD = 9.1) (Fig. S3).
3.1.2 BH14
Two images taken at two distinct locations in the same vein of this sample show a large amount of U (mean of 4 to 11 ppm depending on the image) (Fig. 3a, b) and a very good spread of the pixel U ∕ Pb and Pb ∕ Pb ratios along a linear trend (∼ 2 to 96 for ^{238}U ∕ ^{206}Pb; Fig. 3c, d). Such variation in isotope ratios (and U concentrations) is clearly related to crystal growth zoning (Fig. 3a, b). This composition allows precise dating, as shown by Beaudoin et al. (2018). Indeed, splitting the isotopic maps in pseudospots allows for the calculation of ages of 61.3 ± 3.0 Ma (Fig. 3e) and 61.7 ± 3.0 Ma (Fig. S3), with a MSWD of 0.1 and 0.4 respectively, which are values that are fully consistent with the ages obtained with the robust regressions from the entire maps (Table 1). The latter vary between 60.7 ± 2.0 Ma and 61.8 ± 2.2 Ma according to the diagrams considered (Fig. 3c, d; only TW concordia plots shown). They are identical within the limits of uncertainty to the values obtained with LAICPMS spot analysis, both by Beaudoin et al. (2018) and in our laboratory (63.2 ± 2.2 and 61.2 ± 2.9 Ma respectively; Fig. 3f). They are also similar to the ages calculated by using the approach of Drost et al. (2018) (60.8 ± 2.1 to 61.4 ± 2.1 Ma according to the diagrams considered; Fig. 3g and S3). The age uncertainties calculated by the robust regression are generally less than 1 Ma (without considering systematic uncertainties), which is similar to values obtained by spot analysis and with the method of Drost et al. (2018). This high precision can be related to the very good visual and statistical parameters calculated for this regression (good point alignment, an RSE < 7 % of the Pb_{0} value and a MSWDd < 1.7; Fig. 3c, d, e, f). The colocalisation study shows that the points furthest from the trend line of the robust regression are not located in a same area on the isotopic map. Any selection of a part of the map to possibly improve the age calculations is not justified.
3.1.3 PXG201
This sample is characterised by a moderately high mean U content of 2.5 ppm (Fig. 4a) and by a large spread of the pixel U ∕ Pb and Pb ∕ Pb ratios along a linear trend (∼ 0.5 to 93 for ^{238}U ∕ ^{206}Pb) (Fig. 4b). This favourable composition explains that the ages obtained by the robust regression are comparable within the uncertainties between the TW concordia plot (40.5 ± 1.6 Ma; Fig. 4b), the 86TW plot (38.2 ± 1.6 Ma; Fig. S3) and the isochron plot (40.2 ± 1.6 Ma; Fig. S3) (Table 1). They are also consistent with the age obtained by LAICPMS spot analysis (40.5 ± 2.3 Ma), both in terms of absolute value and of precision (Fig. 4d), and with the age calculated from the pseudospots (40.5 ± 3.3 Ma; Fig. 4c). The slight variations in the central age depending on the diagram used (38.2 to 40.5 Ma) can be related to the statistical parameters of the regression, which are not as good as those obtained for BH14 (MSWDd values from 1.1 to 3, and RSE values between 7 % and 9 % of the Pb_{0} value). This possibly shows some heterogeneity in the isotopic composition – for example, inclusions of detrital material. The colocalisation study shows a globally homogeneous spatial distribution of the values most distant from the regression line in a TW concordia plot. Any inclusions or cement alteration would certainly be better highlighted using a thresholding approach based on trace element contents, as proposed by Drost et al. (2018). The ages calculated with the method of Drost et al. (2018) for the entire map are within uncertainty of the previous values: 42.5 ± 1.8 Ma for the TW concordia plot (Fig. 4e), 40.6 ± 1.7 Ma for the 86TW plot (not shown) and 40.7 ± 1.1 Ma for the isochron plot (not shown). Finally, the ages obtained with both the isotope mapping approaches and LAICPMS spot analysis are centred around ∼ 40 Ma. Such ages are similar to that suggested by the paleomagnetic study of Oms et al. (2003), which pointed to a Bartonian age (∼ 40.0–38.5 Ma). However, they are slightly older (including the limits of uncertainty) than the age proposed by Canudo and Molina (1988) on the basis of the presence of Globigerinatheka semiinvoluta in the sediments containing the sample. This microfauna rather indicates an age between 38.5 and 36.0 Ma according to PremoliSilva et al. (2016) and Wade et al. (2011) and considering the Geologic Time Scale 2012 (GTS2012; Gradstein et al., 2012). The origin of this discrepancy is not yet identified. Nevertheless, the ages obtained are consistent, along with the other proxies, with the transition between marine and continental deposits in the studied area being Bartonian–Priabonian in age.
3.1.4 PXG322
This lake carbonate sample has a very high U content (mean of 17.5 ppm) (Fig. 5a). However, the distribution of the pixel U ∕ Pb and Pb ∕ Pb ratios is more restricted than for the PXG201 sample (Fig. 5b). Accordingly, the age calculated from the pseudospots is less precise than for the previous sample (32.5 ± 3.7 Ma) (Fig. 5c). Conversely, regarding the robust regression, the ages obtained are all comparable and of greater precision (∼ 35.0 ± 2.0 Ma, MSWDd < 1.3 and RSE values lower than 5 % of the Pb_{0} value) (Fig. 5b and S3; Table 1). The ages obtained with the method of Drost et al. (2018) are slightly variable, although they are similar within the limits of the uncertainties (34.2 ± 1.0 Ma to 36.6 ± 1.6, MSWD ∼ 1; only the TW concordia plot is shown in Fig. 5e). Nevertheless, all of the previous ages are in agreement with the age obtained by LAICPMS spot analysis (37.0 ± 3.0 Ma; Fig. 5d) and by the paleomagnetic study of Mochales et al. (2016). The age obtained with the robust regression approach is of higher precision than the age found with LAICPMS spot analysis. However, the age calculated with spot analysis is associated with unsatisfactory statistics (MSWD = 6.3) and could, therefore, be considered as an errorchron. Such a high MSWD is clearly linked with sample heterogeneity as seen from the scatter of error ellipses around the regression line (Fig. 5d). Compositional heterogeneities are also clearly visible with U and Pb elemental maps (Fig. 5a). This is consistent with the heterogenous nature of the sample, composed of pellets and microbial oncoids (Fig. S2). The colocalisation study indicates that there are no obvious links between such heterogeneities and the location of outlier pixel values, which are evenly distributed on the maps and do not justify any subselection of the isotopic map. The age obtained confirms the Bartonian age of the first continental deposits in the eastern part of the Jaca Basin in the southern Pyrenees, bringing an absolute constraint on the end of the nondeposition period above the Boltaña anticline.
3.1.5 JS4 (run 1)
This Carboniferous–Permian lacustrine limestone sample, characterised by high U concentrations (mean ∼ 10 ppm), was previously dated by Drost et al. (2018). These authors obtained precise ages of at 300.8 ± 3.7, 299.6 ± 4.0 and 298.9 ± 5.6 Ma for the TW, 86TW and isochron diagrams respectively. The ages calculated with the robust regression from pixel U ∕ Pb and Pb ∕ Pb ratios (Kerstin Drost, personal communication, 2020) are similar to the previous ages (298.9 ± 3.4, 297.0 ± 3.4 and 300.8 ± 2.8 Ma for the TW, 86TW and isochron diagrams respectively; only the TW concordia plot is shown in Fig. 6a) (Table 1). The very good statistics of the regressions (MSWDd ≤ 1.0, RSE < 5.5 % of the Pb_{0} value) is explained by the large spread of the ratios, defining a clear linear trend in the different plots. Note that to concur with the results given by Drost et al. (2018), these ages are without systematic uncertainties. Adding them leads to final age uncertainties above 8.0 Ma (Fig. 6a).
3.1.6 BM18
Unlike JS4, the BM18 sample, a tectonic calcite vein, is characterised by low U and Pb concentrations (mean of ∼ 200 ppb and ∼ 7 ppb respectively). It was dated precisely to 59.5 ± 2.7 Ma (MSWD = 2.5) by LAICPMS spot analysis, owing to a large, linear spread of the individual ellipses in the TW space (Beaudoin et al., 2018). A similar age (61.0 ± 1.7 Ma in the TW space, MSWD = 1.1) was obtained by Roberts et al. (2020) with the pixelpooling approach of Drost et al. (2018), despite a large scatter of the pixel U ∕ Pb and Pb ∕ Pb isotopic ratios (Fig. 6b). Applying robust regression to these data does not allow the expected age to be reproduced. When systematic uncertainties similar to those of Beaudoin et al. (2018) are considered, the ages obtained are 46.9 ± 2.8 and 48.4 ± 2.5 Ma for the TW (Fig. 6b) and 86TW diagrams (not shown) respectively, which is more than 10 Ma younger than the expected age (Table 1). In contrast, for the isochron diagram, the age of 60.9 ± 2.1 Ma (not shown) is in line with expectations. Added to high RSE values (> 10 % of the Pb_{0} value), as well as MSWDd values ranging from 2.5 to 4.9, these variable ages would lead us to consider the ages as unreliable. The discrepancy between the ages is due to the presence of pixels with high U ∕ Pb ratios, acting as leverage points for the robust regression in the TW and 86TW diagrams.
In order to get closer to the expected age, we added weights to the pixel values used in the robust regression, as is possible with the lmrob library. The weights are based on the density of pixels in each plot, as estimated by a Gaussian kernel density estimate (KDE), the bandwidth of which is estimated by Scott's rule (Scott, 1992). In the case of a large dispersion of the data (in particular linked to the counting statistics, even for a sample of homogeneous age), it is expected that the majority of the points will be grouped along the isochron (i.e. higher density). The use of this approach for the TW and 86TW diagrams gives values similar in uncertainties to that of Beaudoin et al. (2018) but still centred on younger ages (54.4 ± 2.6 and 55.0 ± 2.5 Ma for the TW and 86TW diagrams respectively; Fig. 6c). However, the ages obtained with the three diagrams are not similar within uncertainties, which would again lead to the rejection of the analysis or, as the ages are close to each other, to the consideration that the real age is probably between ∼ 52 and ∼ 63 Ma. Using the approach of Drost et al. (2018), an age of 60.0 ± 2.8 Ma (MSWD = 0.95) is obtained for the TW diagram (Fig. S3), i.e. slightly younger than the value presented in Roberts et al. (2020). Doing the same calculation for the TW86 and isochron diagrams (not presented in Roberts et al., 2020) gives values of 55.6 ± 2.9 Ma (MSWD = 1.0) and 60.1 ± 6.4 Ma (MSWD = 0.5) respectively (not shown). The three ages are centred between 55 and 60 Ma; therefore, they are identical to those obtained with the robust regression, but they are also similar to each other within uncertainties. Thus, for this sample of low U and Pb concentrations, the ages calculated with the approach of Drost et al. (2018) are more reliable than with the robust regression (weighted or not), despite a Pb_{0} value (∼ 0.700) higher than that expected on the basis of spot analyses (∼ 0.590; Beaudoin et al., 2018).
3.2 Unknown samples
3.2.1 ETC2
The isotopic map of this sample (Fig. 7a), taken from a tectonic breccia in the Basque Country (France), was made through several calcite cements as clearly shown in the cathodoluminescence (CL) images (Fig. 7b). The latter are characterised by variable U contents (below the detection limit to > 200 ppm; Fig. 7a). Most pixels have a low U ∕ Pb ratio (< 10), and some dispersion in the ^{207}Pb ∕ ^{206}Pb values for the weakly radiogenic pixels can be noticed (Fig. 7c). However, the U ∕ Pb and Pb ∕ Pb ratios of several pixels define a clear linear trend (Fig. 7c), allowing age calculation. The age obtained from the pseudospots is 100.2 ± 10.1 Ma (Fig. 7d, Table 1). Those obtained with the robust regression are identical, within the limits of uncertainty: 96.1 ± 3.9 Ma and 96.5 ± 4.0 for the respective TW (Fig. 7c) and 86TW plots (Fig. S3) and an age of 95.6 ± 3.3 Ma for the isochron plot (Fig. S3). These similar ages are consistent with the good statistical parameters obtained from the robust regression (MSWDd < 1.2, RSE lower than 8.5 % of the Pb_{0} value). However, both the maps of pixel concentration values and the colocalisation analysis of image data demonstrate that the overdispersion of ^{207}Pb ∕ ^{206}Pb values are located in a restricted area of the map corresponding to one cement (cement C1 on the CL image, Fig. 7b), whereas the highest ^{238}U ∕ ^{206}Pb (i.e. most radiogenic values) that drive the age are located in another restricted area corresponding to another cement (cement C2 on the CL image) (Fig. 7e). Lead by the colocalisation study, the selection of a small area of the map corresponding to cement C2 results in similar but more accurate age calculations than considering the entire map (97.4 ± 3.5 Ma for the TW plot, Fig. 7f; 96.9 ± 3.5 Ma for the 86TW plot, Fig. S3; and 95.5 ± 3.3 Ma for the isochron plot, Fig. S3). The statistical parameters are also better, despite a number of pixels divided by 3 (MSWDd < 0.9, RSE lower than 5 % of the Pb_{0} value). This demonstrates the relevance of the colocalisation approach for age calculation with the robust regression approach, especially when the analytical conditions do not allow trace elements to be analysed at the same time as U and Pb. In contrast, regarding the approach of Drost et al. (2018), the ages obtained from the selected area of the map are less precise than those calculated from the entire map (96.4 ± 6.5 versus 97.9 ± 4.7 Ma for the TW plot, Figs. 7g and S3; 98.0 ± 6.1 versus 100.0 ± 5.0 Ma for the TW86 plot, not shown; and 95.0 ± 7.5 versus 97.1 ± 5.6 Ma for the TW plot, not shown). This is explained by the lower number of ellipses governing the regression.
Finally, in the case of the isotopic image made on the ETC2 sample, the robust regression approach appears to be the most precise. The age obtained indicates cement precipitation during the rifting event that affected the area before the Iberian–Eurasian convergence, which led to the formation of the Pyrenees from the Campanian (Mouthereau et al., 2014).
3.2.2 ARB
ARB sample is a single, homogenous calcite cement partially filling a tectonic breccia. It has a mean U content of ∼ 1.9 ppm (Fig. 8a) and a large spread of the pixel U ∕ Pb and Pb ∕ Pb ratios (∼ 0.2 to ∼ 106.0 for ^{238}U ∕ ^{206}Pb), defining a linear trend (Fig. 8b). Using the robust regression, the calculated ages for this sample are 106.5 ± 3.8 Ma for the TW concordia plot (Fig. 8b), 107.4 ± 3.8 Ma for the 86TW plot (Fig. S3) and 109.7 ± 3.7 Ma for the isochron if the entire image is considered (Table 1; Fig. S3). Again, these values are similar within uncertainties and similar to the age obtained from the pseudospots (108.4 ± 8.1 Ma; Fig. 8d). The slight differences in age between the types of plots (< 3 %) may be related to the high RSE (9 % to 18 % of the Pb_{0} value), which indicates some dispersion around the regression line, despite MSWDd values lower than two. Likewise, the ages obtained by the method of Drost et al. (2018) are slightly variable, although consistent with previous ages within the limits of uncertainties (111.0 ± 4.2 Ma for the TW plot, Fig. 8e; 110.6 ± 3.9 for the 86TW plot, not shown; and 105.3 ± 4.6 Ma for the isochron plot, not shown; with MSWDs of 2.6, 2.5 and 1.1 respectively). For the robust regression, the colocalisation study indicates that among the more distant points from the regression line in the concordia TW plot, almost all correspond to pixels located in the left quarter of the isotope map, where the U and Pb concentrations are the lowest (Fig. 8a). The age uncertainties obtained without considering these pixels are slightly better than those obtained with the entire map: 107.0 ± 3.8 Ma for the TW concordia plot (Fig. 8c), 107.8 ± 3.8 Ma for the 86TW plot (Fig. S3) and 109.7 ± 3.7 Ma for the isochron plot (Fig. S3). However, contrary to ETC2, the statistical parameters are not better (MSWDd between 1.4 and 2.4, and RSE between 9.5 and 18.0 % of the Pb_{0} value) due to the disappearance of several pixels located close to the robust regression line. In contrast, applying the method of Drost et al. (2018) to the map subset results in better age uncertainties and MSWD values compared with the entire isotopic map, with values of 110.5 ± 4.0 Ma for the TW concordia plot (MSWD = 2.0; Fig. S3), 110.4 ± 3.8 Ma for the 86TW plot (MSWD = 1.1; not shown) and 104.7 ± 4.4 Ma for the isochron plot (MSWD = 0.6; not shown).
Taken together, the ETC2 and ARB samples suggest that the emplacement of tectonic breccia in the western Pyrenees was related to extensive synrift fault development, rather than to the subsequent tectonic inversion.
3.2.3 C6265D5
This sample is a dolomite vein from the northern Pyrenees, presumably formed in a tectonic context identical to ETC2 and ARB samples. Its mean U and Pb contents are low (230 and 60 ppb respectively; Fig. 9a), and a large scatter of the U ∕ Pb and Pb ∕ Pb pixel isotope ratios can be evidenced (Fig. 9b). As a result, the age calculated from the pseudospots has large uncertainties (115.2 ± 17.4 Ma; Fig. 9c; Table 1), and the use of robust regression does not make it possible to obtain consistent ages between the three types of plots (84.5 ± 4.7, 82.4 ± 4.4 and 105.3 ± 3.6 Ma for the TW, 86TW and isochron plots respectively; Figs. 8b and S3). These ages are associated with MSWDd values close to unity but very high RSE values (21 % to 28 % of the Pb_{0} value). Weighting the regression by KDE for the three plots (according to the method presented for sample BM18) allows one to resolve these inconsistencies (97.0 ± 5.7, 94.5 ± 5.3 and 93.7 ± 4.5 Ma for the TW, 86TW and isochron plots respectively) and to get better statistical parameters (MSWDd = 0.8 to 1, and RSE = 10 % to 21 % of the Pb_{0} value) (Figs. 9d, S3). However, the new ages are outside the uncertainty of those obtained from the pseudospots (i.e. younger), which calls their reliability into question. On the other hand, the ages calculated by the method of Drost et al. (2018) are also comparable to those calculated from the pseudospots, with values of 110.3 ± 6.2, 102.7 ± 6.2 and 104.7 ± 7.8 for the TW (Fig. 8e), 86TW and isochron plots respectively and with MSWD values lower than 1.8. It therefore seems that, in a manner analogous to BM18, the results obtained by robust regression for C6265D5 are less reliable than those obtained by the pixelpooling approach of Drost et al. (2018), except for the isochron plots. Nevertheless, although they are obtained on dolomite and therefore possibly biased by a matrix effect (Elisha et al., 2020), these ages are consistent with those obtained for samples ETC2 and ARB, which once again suggests a link with the rifting event recorded in the Pyrenees during the Cretaceous.
4.1 Interest of dating from isotope maps
The pioneering work of Drost et al. (2018) has shown the obvious value of the isotope map dating approach, where data on minor and trace element concentrations are used to select the most suitable pixels for age determination. This initial step seems essential and should be followed up systematically if the analytical conditions allow for it (e.g. use of an ICPMS with a fast mass scan such as a quadrupole or time of flight). If not, the operator must remain attentive to the quality of the data by checking the scattering of the pixels around the regression lines and the possible presence of phase mixtures. The use of a colocalisation approach as presented in our study can be an alternative way, although less efficient, of carrying out this careful study of the spatial distribution of the furthest pixel values from the regression line. However, the examples presented above, in particular that of the ETC2 sample, show that colocalisation is more suited to the robust regression approach than to that of pixel pooling. In any case, it is fundamental to carry out an initial petrographic characterisation of the samples as recalled by Roberts et al. (2020). Taking these precautions into account, our results and those of Drost et al. (2018) and Roberts et al. (2020) show that the ages obtained from the imagemapping approach are perfectly comparable, both in value and uncertainty, with those obtained from more common methods (LAICPMS spot analysis and isotope dilution). Additionally, it should be pointed out that this approach has several advantages in comparison with the more common ones. First, the imagemapping approach is particularly suitable for easily dating multiple generations of cements filling porosities or fractures. Such an approach is fast and simple to set up, as it is sufficient to map a large portion of the sample while covering several generations of cements with a number of pixels high enough for age calculation; the user can then select the pixels (preferentially by a geochemical approach) corresponding to the different cements to calculate the ages. This methodology has been carried out successfully on sample ARB with colocalisation. Second, contrary to spot analyses, the simultaneous mapping of trace elements useful for both geochemical characterisation and for dating makes it possible to gather all of the necessary information in a single analysis, without changing the analytical parameters of the LAICPMS, which may save a considerable amount of time. Third, a complete isotope map can be acquired very quickly (less than 30 min, without the reference materials). This allows one to consider the analysis of many samples over a limited number of analytical sessions. Therefore, it is likely that this method of U–Pb dating of carbonates will develop rapidly in the geochronology community, alongside the now wellestablished LAICPMS spot analysis approach.
4.2 Benefits and limitations of the robust regression to U–Pb image dating
The data processing method presented here, namely robust regression through the isotope ratio values of the pixels represented in TW or isochron diagrams, is simple and allows a fast age determination including its uncertainty. We recall, as previously mentioned, the importance of a preliminary selection of these pixels on the basis of the concentrations of trace elements of the carbonates analysed. As an alternative to the datapooling approach presented by Drost et al. (2018), it has the advantage of not involving the operator in the choice of the number of pixels in each pseudoanalysis, as all pixels are considered for a single regression. Their representation in the appropriate diagrams gives a quick idea of the quality of the data, assisted by statistical parameters such as the RSE or MSWDd. Note that the interest in using robust regression for U–Pb dating was already underlined by Powell et al. (2002) and very recently by Powell et al. (2020). It is also implemented in the famous Isoplot software (Ludwig, 2012). Ludwig (2012) stressed the precautions to be taken when using robust regression to construct isochrons, as the absence of uncertainties on each point makes it impossible to distinguish the share of scatter around the regression due to “geological/geochemical complications” from the scatter due to analytical errors alone. In conventional spot analyses, the MSWD calculation is an efficient way to determine the share of each source of uncertainty. Its use in the context of dating from isotopic images, whether with the approach of Drost et al. (2018) or our own, does not solve this problem as the MSWD is calculated from an “artificial” grouping of pixels whose associated ellipses do not represent analytical uncertainties. At most, the MSWD (or MSWDd) is here a means of verifying the correct alignment of values in the TW and isochron diagrams, with values close to one indicating good alignment, information that can be complemented by the moving average calculation as we propose, among other possibilities. For example, the very recent approach proposed by Powell et al. (2020), with calculation of a normalised (robust) median of the absolute values of the residuals, could provide a suitable means to distinguish isochrons from errorchrons. This still has to be implemented in our data reduction scheme.
Among the negative points of the robust regression, the example of samples BM18 and C6265D5 shows the sensitivity of the method to the dispersion of the U ∕ Pb and Pb ∕ Pb ratios. For this type of lowU sample, where noisy data are linked to low counting statistics, the ages calculated without weighting are too young for the TW and 86TW plots. They are also inconsistent between the TW, 86TW and isochron plots. This would lead to rejecting samples which, with the method proposed by Drost et al. (2018), are validated by obtaining consistent ages not only between the different diagrams but also with those obtained by spot or by pseudospot analyses. Weighting the robust regression by KDE improves the results by giving consistent ages between the diagrams. However, the comparison with the expected ages shows that the ages calculated in this way are too young by around 5 % to 10 % (55–60 Ma instead of 60 Ma for BM18, and 94–97 Ma instead of 118 Ma for C6265D5). For such samples, the tests carried out show the greatest reliability of the method of Drost et al. (2018) for imagebased dating. However, as this latter approach might give erroneous Pb_{0} values (example of BM18), we argue that LAICPMS spot analyses should be preferred when U concentrations are far below parts per million.
U–Pb LAICPMS carbonate dating from isotope images is an extremely recent approach that, compared with more common dating methods such as LAICPMS spot analyses, has the great advantage of allowing the use of images for both sample screening and dating. In this contribution, U–Pb ages of several carbonate samples were obtained from micronscale LAICPMS isotopic maps, using a new, alternative data processing method based on robust regression. The method, which is very simple to set up, consists of extracting the U ∕ Pb and Pb ∕ Pb pixel values from the image, plotting them in a diagram used for dating (such as a concordia diagram) and performing robust regression to calculate the age. Prior to dating, we also conducted a colocalisation study of the sample isotope ratios that was used to locate the most suitable areas for age determination. For most samples that had already been dated by U–Pb geochronology using conventional LAICPMS spot analysis or isotope dilution, the robust regression approach method provided similar ages and sometimes better uncertainties compared with both others methods, with an ablation duration of less than 40 min. For lowU samples with noisy U ∕ Pb and Pb ∕ Pb pixel values, however, the ages were biased towards overly young values by 5 % to 10 %. In this case, the existing imagebased dating approach of Drost et al. (2018) gave better results. Thus, the robust regression processing of LAICPMS isotope map data appears to be a simple means to proceed to carbonate U–Pb geochronology, provided that a careful examination of the pixel ratio data is performed beforehand.
The methodology for LAICPMS spot analyses, the screening maps and the plots used for age calculations are presented in the Supplement. The pixel values of the isotopic images for the masses 238, 232, 208, 207 and 206 are presented in Table S1 in the Supplement.
The supplement related to this article is available online at: https://doi.org/10.5194/gchron3672021supplement.
GH, FC and CPe designed the experiments. GH, FC, CPe, PAG and GM carried them out. GH and CPa developed the model code and performed the simulations. OC and JPG supervised the project and provided the financial support for most experiments. GH prepared the paper with contributions from all coauthors.
The authors declare that they have no conflict of interest.
This article is part of the special issue “In situ carbonate U–Pb geochronology”. It is a result of the Goldschmidt conference, Barcelona, Spain, 18–23 August 2019.
The authors are grateful to Gaëlle Barbotin for help during the analytical sessions. Jose Ignacio Canudo is warmly thanked for a discussion on the comparison of U–Pb and biostratigraphic ages in the Jaca Basin. The initial paper was greatly improved by comments from the two anonymous reviewers.
This work benefited from funding from the Universite de Pau et des Pays de l'Adour (AAP Incitatif Recherche) and from Total SA (project DAC).
This paper was edited by Axel Gerdes and reviewed by two anonymous referees.
Aramendía, M., Rello, L., Bérail, S., Donard, A., Pécheyran, C., and Resano, M.: Direct analysis of dried blood spots by femtosecondlaser ablationinductively coupled plasmamass spectrometry. Feasibility of splitflow laser ablation for simultaneous trace element and isotopic analysis, J. Anal. Atom. Spectrom., 30, 296–309, 2015.
Ballihaut, G., Claverie, F., Pécheyran, C., Mounicou, S., Grimaud, R., and Lobinski, R.: Sensitive Detection of Selenoproteins in Gel Electrophoresis by High Repetition Rate Femtosecond Laser AblationInductively Coupled Plasma Mass Spectrometry, Anal. Chem., 79, 6874–6880, 2007.
Beaudoin, N., Lacombe, O., Roberts, N. M. W., and Koehn, D.: UPb dating of calcite veins reveals complex stress evolution and thrust sequence in the Bighorn Basin, Wyoming, USA, Geology, 46, 1015–1018, https://doi.org/10.1130/G45379.1, 2018.
Canudo, J. I. and Molina, E.: Biocronología de foraminíferos planctónicos de la secuencia deposicional de Jaca (Pirineo Aragonés): Eoceno medio y superior, in: Proceedings of the II Congreso Geológico de España, Granada, Spain, 273–276, 1988 (in Spanish).
Drost, K., Chew, D., Petrus, J. A., Scholze, F., Woodhead, J. D., Schneider, J. W., and Harper, D. A. T: An image mapping approach to UPb LAICPMS carbonate dating and applications to direct dating of carbonate sedimentation, Geochem., Geophy., Geosy., 19, 46314648, https://doi.org/10.1029/2018GC007850, 2018.
Elisha, B., Nuriel, P., KylanderClark, A., and Weinberger, R.: Towards insitu U–Pb dating of dolomites, Geochronology Discuss. [preprint], https://doi.org/10.5194/gchron202019, in review, 2020.
Godeau, N., Deschamps, P., Guihou, A., Leonide, P., Tendil, A., Gerdes, A., Hamelin, B., and Girard, J.P.: UPb dating of calcite cement and diagenetic history in microporous carbonate reservoirs: Case of the Urgonian Limestone, France, Geology, 46, 247–250, https://doi.org/10.1130/G39905.1, 2018.
Gradstein, F. M, Ogg, J. G., Schmitz, M. D., Ogg, G. M., Agterberg, F. P., Anthonissen, D. E., Becker, T. R., Catt, J. A., Cooper, R. A., Davydov, V. I., Gradstein, S. R., Henderson, C. M., Hilgen, F. J., Hinnov, L. A., McArthur, J. M., Melchin, M. J., Narbonne, G. M., Paytan, A., Peng, S., PeuckerEhrenbrink, B., Pillans, B., Saltzman, M. R., Simmons, M. D., Shields, G. A., Tanaka, K. L.,Vandenberghe, N., Van Kranendonk, M. J., Zalasiewicz, J., Altermann, W., Babcock, L. E., Beard, B. L., Beu, A. G., Boyes, A. F., Cramer, B. D., Crutzen, P. J., van Dam, J. A., Gehling, J. G., Gibbard, P. L., Gray, E. T., Hammer, O., Hartmann, W. K., Hill, A. C., Paul F. Hoffman, P. F., Hollis, C. J., Hooker, J. J., Howarth, R. J., Huang, C., Johnson, C. M., Kasting, J. F., Kerp, H., Korn, D., Krijgsman, W., Lourens, L. J., MacGabhann, B. A., Maslin, M. A., Melezhik, V. A., Nutman, A. P., Papineau, D., Piller, W. E., Pirajno, F., Ravizza, G. E., Sadler, P. M., Speijer, R. P., Steffen, W., Thomas, E., Wardlaw, B. R., Wilson, D. S., and Xiao, S.: he Geologic Time Scale 2012, Elsevier, Boston, USA, 1129 pp., 2012.
Guillong, M., Wotzlaw, J.F., Looser, N., and Laurent, O.: Evaluating the reliability of U–Pb laser ablation inductively coupled plasma mass spectrometry (LAICPMS) carbonate geochronology: matrix issues and a potential calcite validation reference material, Geochronology, 2, 155–167, https://doi.org/10.5194/gchron21552020, 2020.
Hill, C. A., Polyak, V. J., Asmerom, Y., and Provencio, P. P.: Constraints on a Late Cretaceous uplift, denudation, and incision of the Grand Canyon region, southwestern Colorado Plateau, USA, from UPb dating of lacustrine limestone, Tectonics, 35, 896–906, https://doi.org/10.1002/2016TC004166, 2016.
Horstwood, M. S. A., Košler, J., Gehrels, G., Jackson, S. E., Mclean, N. M., Paton, C., Pearson, N. J., Sircombe, K., Sylvester, P., Vermeesch, P., Bowring, J. F., Condon, D. J., and Schoene, B.: CommunityDerived Standards for LAICPMS U(Th)Pb Geochronology – Uncertainty Propagation, Age Interpretation and Data Reporting, Geostand. Geoanal. Res., 40, 311332, https://doi.org/10.1111/j.1751908X.2016.00379.x, 2016.
Koller, M. and Stahel, W. A.: Sharpening Waldtype inference in robust regression for small samples, Comput. Stat. Data An., 55, 2504–2515, 2011.
Koller, M. and Stahel, W. A.: Nonsingular subsampling for regression Sestimators with categorical predictors, Computation. Stat., 32, 631–646, https://doi.org/10.1007/s001800160679x, 2017.
Li, Q., Parrish, R. R., Horstwood, M. S. A., and McArthur, J. M.: UPb dating of cements in Mesozoic ammonites, Chem. Geol., 376, 76–83, 2014.
Ludwig, K. R.: User's manual for Isoplot version 3.75–4.15: a geochronological toolkit for Microsoft Excel, Berkeley Geochronological Center Special Publication, Berkeley, CA, USA, 2012.
Maechler, M., Rousseeuw, P., Croux, C., Todorov, V., Ruckstuhl, A., SalibianBarrera, M., Verbeke, T., Koller, M., Conceicao, E. L. T., and Anna di Palma, M.: robustbase: Basic Robust Statistics R package version 0.935, available at: http://CRAN.Rproject.org/package=robustbase (last access: 1 April 2020), 2019.
Mochales, T., Pueyo, E. L., Casas, A. M., and Barnolas, A.: Restoring paleomagnetic data in complex superposed folding settings: the Boltaña anticline (Southern Pyrenees), Tectonophysics, 671, 281–298, 2016.
Mouthereau, F., Filleaudeau, P. Y., Vacherat, A., Pik, R., Lacombe, O., Fellin, M. G., Castelltort, S., Christophoul, F., and Masini, E.: Placing limits to shortening evolution in the Pyrenees: role of margin architecture and implications for the Iberia/Europe convergence, Tectonics, 33, 2283–2314, https://doi.org/10.1002/2014TC003663, 2014.
Nuriel, P., Craddock, J., KylanderClark, A. R., Uysal, T., Karabacak, V., Dirik, R. K., Hacker, B. R., and Weinberger, R.: Reactivation history of the North Anatolian fault zone based on calcite agestrain analyses, Geology, 47, 465–469, 2019.
Oms, O., DinarèsTurell, J., and Remacha, E.: Magnetic stratigraphy from deep clastic turbidites: An example from the Eocene Hecho Group (southern Pyrenees), Studia Geophysica et Geodetica, 47, 275–288, 2003.
Parrish, R. R., Parrish, C. M., and Lasalle, S.: Vein calcite dating reveals Pyrenean orogen as cause of Paleogene deformation in southern England, J. Geol. Soc. London, 175, 425–442, https://doi.org/10.1144/jgs2017107, 2018.
Powell, R., Woodhead, J., and Hergt, J.: Improving isochron calculations with robust statistics and the bootstrap, Chem. Geol., 185, 91–204, 2002.
Powell, R., Green, E. C. R., Marillo Sialer, E., and Woodhead, J.: Robust isochron calculation, Geochronology, 2, 325–342, https://doi.org/10.5194/gchron23252020, 2020.
Premoli Silva, I., Wade, B. S., and Pearson, P. N.: Taxonomy, biostratigraphy, and phylogeny of Globigerinatheka and Orbulinoides, in: Atlas of Eocene Planktonic Foraminifera, edited by: Pearson, P. N., Olsson, R. K., Hemleben, C., Huber, B. T., and Berggren, W. A., Cushman Foundation for Foraminiferal Research, Fredericksburg, USA, Special Publication, 41, 169–212, 2016.
Ring, U. and Gerdes, A.: Kinematics of the AlpenrheinBodensee graben system in the Central Alps: Oligocene/Miocene transtension due to formation of the Western Alps arc, Tectonics, 35, 1367–1391, https://doi.org/10.1002/2015TC004085, 2016.
Roberts, N. M. W. and Walker, R. J.: UPb geochronology of calcitemineralized faults: Absolute timing of riftrelated fault events on the northeast Atlantic margin, Geology, 44, 531–534, 2016.
Roberts, N. M. W., Rasbury, E. T., Parrish, R. R., Smith, C. J., Horstwood, M. S. A., and Condon, D. J.: A calcite reference material for LAICPMS UPb geochronology, Geochem., Geophy., Geosy., 18, 2807–2814, https://doi.org/10.1002/2016GC006784, 2017.
Roberts, N. M. W., Drost, K., Horstwood, M. S. A., Condon, D. J., Chew, D., Drake, H., Milodowski, A. E., McLean, N. M., Smye, A. J., Walker, R. J., Haslam, R., Hodson, K., Imber, J., Beaudoin, N., and Lee, J. K.: Laser ablation inductively coupled plasma mass spectrometry (LAICPMS) U–Pb carbonate geochronology: strategies, progress, and limitations, Geochronology, 2, 33–61, https://doi.org/10.5194/gchron2332020, 2020.
Scott, D. W.: Multivariate Density Estimation: Theory, Practice, and Visualization, John Wiley and Sons, New York and Chicester, USA, 1992.
Vermeesch, P.: IsoplotR: A free and open toolbox for geochronology, Geosci. Front., 9, 1479–1493, https://doi.org/10.1016/j.gsf.2018.04.001, 2018.
Wade, B. S., Pearson, P. N., Berggren, W. A., and Pälike, H.: Review and revision of Cenozoic tropical planktonic foraminiferal biostratigraphy and calibration to the geomagnetic polarity and astronomical time scale, EarthSci. Rev., 104, 111–142, 2011.
Woodhead, J. D. and Hergt, J. M.: Strontium, neodymium and lead isotope analyses of NIST glass certi?ed reference materials: SRM 610, 612, 614, Geostandard. Newslett., 25, 261–266, https://doi.org/10.1111/j.1751908X.2001.tb00601.x, 2001.
Zeitvogel, F. and Obst, M.: ScatterJn: An ImageJ Plugin for ScatterplotMatrix Analysis and Classification of Spatially Resolved Analytical Microscopy Data, J. Open Res. Softw., 4, e5, https://doi.org/10.5334/jors.89, 2016.