PbIso: an R package and web app for calculating and plotting Pb isotope data

Lead (Pb) isotopes are used in a range of science applications, including plate tectonics, studies of early Earth evolution, and archaeometry. Pb isotopes have been used in geology since the early 1900s (Davis et al. 2003). The U–Th–Pb isotope system is based on the decay of two isotopes of uranium to two different isotopes of Pb and the decay of Th to a third isotope of Pb:

• ²³⁸U → ²⁰⁶Pb   half-life ∼4468 million years.

• ²³⁵U → ²⁰⁷Pb   half-life ∼700 million years.

• ²³²Th → ²⁰⁸Pb   half-life ∼14 000 million years.

A fourth isotope of Pb, ²⁰⁴Pb, which is of primordial origin and not produced by radioactive decay, is used to create isotope ratios (²⁰⁶Pb/²⁰⁴Pb, ²⁰⁷Pb/²⁰⁴Pb, and ²⁰⁸Pb/²⁰⁴Pb), which facilitate measurements using mass spectrometers. A major advantage of this combined U–Th–Pb isotope system is that three decay chains should provide internally consistent results, a feature that is not available for other isotope systems. Decay constants for the three radioactive nuclides are very different, as are Th/U, U/Pb, and Th/Pb ratios in many minerals and rocks, leading to many interesting uses for these isotope systems.

Two principal foci of endeavour based on the U–Th–Pb isotope system have been in (1) geochronology, primarily using minerals with elevated U/Pb ratio, and (2) petrogenesis.

The first of these mostly uses minerals with elevated U/Pb or Th/Pb such as zircon, baddeleyite, titanite, and monazite in which the radioactive decay of uranium or thorium produce large ²⁰⁶Pb/²⁰⁴Pb, ²⁰⁷Pb/²⁰⁴Pb, or ²⁰⁸Pb/²⁰⁴Pb ratios and is therefore generally referred to as the radiogenic U–Th–Pb system. Petrogenetic studies, in contrast, tend to concentrate on either whole rock samples or on minerals with fairly low to negligible U/Pb or Th/Pb, and this version of the U–Th–Pb isotope system is thus typically referred to as the “common Pb” system. The large difference in decay constant, and hence half-life, for ²³⁸U and ²³⁵U produces very different decay characteristics for old (Archean) rocks and minerals when compared with younger samples. Since ²³⁵U decays much more rapidly than ²³⁸U, the ²⁰⁷Pb/²⁰⁴Pb ratio increases much more rapidly than ²⁰⁶Pb/²⁰⁴Pb in early Earth history. At younger (Mesoproterozoic to Phanerozoic) times, there is little change in ²⁰⁷Pb/²⁰⁴Pb, while ²⁰⁶Pb/²⁰⁴Pb continues to increase (Fig. 1).

As with all isotope systems, complexity may be introduced by geological processes such as multiple disturbance events or mixing of radiogenic products. Unlike some of the other petrogenetic isotope systems routinely used today (e.g., Lu–Hf and Sm–Nd) and in the past (Rb–Sr) that exclusively link to silicate petrogenetic processes, the U–Th–Pb isotope system provides insights for both silicate and sulfide processes. Although not frequently used, the Re–Os isotope system also provides information on sulfide petrogenetic processes but is not directly affected by silicate processes. U, Th, and Pb concentrations are greater in crustal material than in most mantle lithologies, and therefore the Pb isotope system is sensitive to interaction (contamination or partial melting) of crust and introduction of this material to mantle-derived mafic lithologies.

The evolution of Pb isotope ratios through time for various conceptual Earth reservoirs is captured by a number of growth models. The simplest assumes a primordial Pb isotope composition equivalent to that of the Canyon Diablo meteorite (Blichert-Toft et al. 2010) and growth of the Pb isotope ratios to values representative of some modern environment. This model was first developed by Holmes (Holmes 1946) and by Houtermans (Houtermans 1946) and is generally referred to as either a single stage model or as the Holmes–Houterman model. It is known that this model does not provide a good representation of Pb isotope evolution for the Earth. More complex two-stage models have been developed, the most frequently used being that of Stacey and Kramers (1975). This model was based on the Pb isotope composition of a number of ore lead minerals and is more accurately described as representing well-mixed ore environments. It had a first stage from the time of formation of the Earth (as represented by Canyon Diablo troilite) until 3.7 Ga, when the U/Pb and Th/Pb composition changed. Models such as that of Cumming and Richards (1975) and Tolstikhin et al. (2006) are more complex, with continually changing U/Pb ratios. Other models are defined for limited geographic regions (e.g., Thorpe (1999) for the Abitibi Belt of Canada) or so recently published that they have not yet been extensively used (e.g. Maltese and Mezger 2020). Each of these two- or multistage models essentially describes a single evolution curve represented by a particular reservoir, although most are generally thought of as representing a version of “average Earth”. In some cases, alternative curves with different U/Pb and Th/Pb composition may be calculated based on the age and composition of the second stage of the original model. Zartman and colleagues (Doe 1979; Zartman and Doe 1981; Zartman and Richardson 2005) developed a model describing the evolution of several distinct reservoirs in the Earth, specifically lower crust, upper crust, mantle, and orogene (a mix of material from the other three reservoirs). This model closely resembles that of Stacey and Kramers (1975) and Cumming and Richards (1975).

Original application of these various models was related to the study of U–Th–Pb in whole rock silicate systems as an additional tracer, supplementing Rb–Sr and Sm–Nd studies. It has, however, been recognised for some time that the various isotope systems may be decoupled from each other, either because of fluid–rock interaction processes (Johnson and DePaolo 1994; Eglington 2019) or because the magmas involved tap material for each isotope system from multiple, distinct reservoirs (Vervoort et al. 1994; Stracke 2012). The viability of the latter has become increasingly important with recent suggestions that distinct zones of concentration of sulfide minerals may occur in the crust and mantle lithosphere that may be tapped by younger mantle melts (Holwell et al. 2022) and the recognition that there may be distinct spatial variations in Pb isotope composition of different ore deposits and crustal domains (Gulson 1986; Warren and Thorpe 1994; Carr et al. 1995; Luais and Hawkesworth 2002; Be’eri-Shlevin et al. 2010; Champion and Huston 2016; Huston et al. 2016; Eglington 2018a; Halla 2018).

Since ²⁰⁶Pb/²⁰⁴Pb or ²⁰⁷Pb/²⁰⁴Pb ratios are different at different times, as is the case for epsilon values of Nd or Hf isotopes, it does not help to plot these values on maps. A more distinctive parameter is the model source ²³⁸U/²⁰⁴Pb (μ) that is a function of model starting composition, model age and event age (mineralisation or crust-forming episode), and the isotope composition of the material produced. The ²³⁸U/²⁰⁴Pb ratio needed to derive an initial composition for one or more samples may be calculated by rearranging the uranogenic isochron equations and solving iteratively if the age of the samples is known (Gale and Mussett 1973; Harmer et al. 1995; Andersen 1998; Albarède et al. 2012; Champion and Huston 2016; Eglington 2018a; Halla 2018). The model source μ (²³⁸U/²⁰⁴Pb) derived in this way may be thought of as the composition of the source region from which rocks or minerals are produced. Calculation of this parameter is quite routine for aficionados of Pb isotopes but is not routinely understood by most geologists, even though its use is again becoming more common in the ore deposit research community. A variety of software have been produced over time for use with isotope data. Some are primarily intended for calculating isochrons, regression ages, and initial compositions (Andersen 1998; Ludwig 2001; Eglington 2018b; Vermeesch 2018). Some have explicit capabilities for model source μ calculations, and the form that these occur in has varied through time as a function of the preferred software development language or environment at the time. Ludwig (2001) and Andersen (1998) provided ways to perform the calculations in Excel spreadsheets, Eglington (2018b) using a stand-alone graphical user package, and Gaab et al. (2006) using the Octave environment. Here, we provide a standalone R package and R code to perform the calculations since R is becoming increasingly popular for scientific investigations.

More comprehensive introductions to Pb isotope plots and their interpretation have been provided elsewhere and are not repeated here (Gale and Mussett 1973; Faure 1977; Halla 2018).

In recent years, several tools have been developed and adapted for various isotopic and geochemical datasets, including, IsoplotR (Vermeesch 2018) (an R implementation of the Excel Isoplot plugin (Ludwig 2001)), provenance (Vermeesch et al. 2016), and detzrcr (Andersen et al. 2018) in R, and pyrolite (Williams et al. 2020) in Python. Some of these packages have a linked graphic user interface, which makes them accessible to users of various programming experience. The power of these tools is the ability to apply them to large datasets and integrate them with other powerful statistical and visualisation packages, which is becoming increasingly important as many disciplines within Earth sciences involve big data analytics.

Typically, Pb isotope calculations and Pb evolution models are performed in makeshift spreadsheets with little transparency of how they are actually calculated. The methods for calculating these different values also vary among publications, often with poor documentation, making reproducibility difficult. Our aim is to provide a review of Pb isotope systematics and how these are incorporated into various Pb isotope calculations. We document the different calculations used for Pb isotope data and Pb evolution models and how they have been implemented into the PbIso R package. The PbIso functions allow for simplicity in only requiring minimal Pb isotope measurements as inputs, while also allowing users the flexibility of setting optional input values or defining their own Pb evolution models.

PbIso is intended for users interested in modelling the evolution of various systems from Pb isotopes, such as calculating model age, model source μ, and initial isotope ratios. These calculations are particularly applicable to a wide range of ore deposit studies (e.g., Huston et al. 2010), and plate tectonic studies (e.g., Blichert-Toft et al. 2016). PbIso also allows rapid modelling of user-defined Pb evolution curves, which is important for understanding Earth evolution as well as the evolution of many ore-forming regions around the world. For individual samples, we recommend using the functions in IsoplotR (Vermeesch 2018), which allows details such as factoring in uncertainties, performing regressions, and sample statistics to be calculated.

The decay of ²³⁸U to ²⁰⁶Pb; ²³⁵U to ²⁰⁷Pb; and ²³²Th to ²⁰⁸Pb, and the non-radiogenic ²⁰⁴Pb can be used to understand the history of certain crustal and mantle processes of a mineral or rock. The abundance of the radiogenic isotopes increases with time (Fig. 2), so samples with significantly different ages cannot easily be compared using the ²⁰⁶Pb/²⁰⁴Pb, ²⁰⁷Pb/²⁰⁴Pb, and ²⁰⁸Pb/²⁰⁴Pb ratios. To allow easy comparison of samples across broad time periods, the model source μ (²³⁸U/²⁰⁴Pb) and κ (²³²Th/²³⁸U) values can be used, calculated using an independently constrained age for each sample. These also allow us to compare samples/deposits of varying ages to modelled reservoirs such as upper crust, mantle, and lower crust (Zartman and Doe 1981); Bulk Silicate Earth (Maltese and Mezger 2020); or region specific models (e.g., Abitibi–Wawa in Superior Province (Thorpe 1999)).

PbIso includes functions for calculating initial Pb isotope ratios, model age, model source μ (²³⁸U/²⁰⁴Pb) and time-integrated κ (²³²Th/²³⁸U), as well as plotting parameters such as model curves, paleoisochron lines and y-intercepts, and isochrons. Using the package within R allows flexibility in applying the functions to the user’s own datasets and the ability to use the wide array of plotting and statistical tools available in R.

We have also implemented some of the basic features of PbIso into an online Shiny R application (see https://shereearmistead.github.io/software/pbiso), which requires no knowledge of R, making it accessible for users. The app allows input of user data, including sample name/ID, age, ²⁰⁶Pb/²⁰⁴Pb, ²⁰⁷Pb/²⁰⁴Pb, and ²⁰⁸Pb/²⁰⁴Pb ratios. Users can then export the processed data as a .xlsx file, which will include the calculated columns for initial Pb isotope ratios, model age, model source μ (²³⁸U/²⁰⁴Pb), and time-integrated κ (²³²Th/²³⁸U), as well as a sheet containing the model parameters. Four basic plots are automatically generated in the Shiny application based on the user data, and can be downloaded as .pdf figures.

The functions in PbIso take one or more of the basic input parameters t (time (Ma)), x (²⁰⁶Pb/²⁰⁴Pb), y (²⁰⁷Pb/²⁰⁴Pb), and z (²⁰⁸Pb/²⁰⁴Pb) to perform the calculations. For advanced usage, the functions can also optionally take the values for different model parameters (summarised in Table 1). The calculations and functions used in PbIso assume a starting composition and model following Stacey and Kramers (1975) second stage model, although this can easily be changed if an alternative model (e.g., Stacey and Kramers (1975) single stage; Maltese and Mezger (2020) Bulk Silicate Earth; or others) is preferred.

The package can be downloaded from https://github.com/ShereeArmistead/PbIso or can be installed within R by running the following:

Note: the first two lines only need to be run once to install the package on a user’s computer. The third line needs to be run every time a user wants to use PbIso in a new R session.

All of the functions in PbIso are designed for ease of use while also allowing flexibility in changing model parameters. The required inputs are outlined in subsequent sections of this manuscript, but we have included a brief overview of the different ways these functions can be used below, using the Calc64() function as an example, which only requires one input (age). The formatted code in the following sections includes the input line of code (e.g., Calc64(2700) in the example below), and the output value given by R (e.g., 13.63662 in the example below), indicated by the line beginning with [1].

The most basic usage is to simply include the one required input parameter, in this case age:

As stated previously, the default starting parameters are based on Stacey and Kramers (1975) second stage model; however, these can be manually overridden by specifying them in the function. The optional parameters, in this case, T1, X1, and Mu1, can be specified as:

Not all of the optional parameters need to be defined. For example, accepting the defaults for T1 and X1, but modifying Mu1 to 8:

Table 2 summarises the PbIso functions and their required and optional input parameters. See Table 1 for the descriptions and default values for these model parameters. Note that the decay constants are included as optional arguments in PbIso functions for maximum flexibility; however, we advise against modifying these unless there are good constraints on alternative values.

Rather than setting the optional parameters manually for every calculation, a predefined model can be used. We have incorporated several commonly used models (see Table 3), for example, the Stacey and Kramers (1975) first stage model.

Or the Maltese and Mezger (2020) Bulk Silicate Earth model:

Alternatively, users may define starting parameters for their own Pb evolution models, which is particularly useful if this is needed for multiple calculations, and to generate model curves (see Section 4). To define your own model, use the list() function in R to define the starting parameters:

Note that for Calc64(), only age, T1, X1, and Mu1 are required (not Y1, Z1, or W1); however, to make the user-defined model applicable to other PbIso functions, it is best to include parameters required for all functions of interest.

Models that we have included as options in PbIso are given in Table 3.

Note that ω (²³²Th/²⁰⁴Pb) in eq. 4 can also be expressed as μ ⋅ κ, which is equivalent to ²³⁸U/²⁰⁴Pb ⋅ ²³²Th/²³⁸U.

The above equations are implemented in PbIso by the functions Calc64(), Calc74(), and Calc84(), respectively. These functions can be used in a number of ways. For simply calculating the value of each isotope ratio on the Stacey and Kramers (1975) average ore lead curve at a given time, only the age is required.

For example, the ²⁰⁶Pb/²⁰⁴Pb, ²⁰⁷Pb/²⁰⁴Pb, and ²⁰⁸Pb/²⁰⁴Pb ratios on the Stacey and Kramers (1975) curve at a given time, say 2700 Ma, is given by

See Table 2 for the optional parameters for each of these functions to allow a different Pb evolution model to be used.

It is not possible to solve this equation directly for the model age (t_sample), so a Newton–Raphson iterative calculation is implemented using the uniroot() function in the stats package in R. This is implemented in PbIso using the CalcModAge() function. Only the ²⁰⁶Pb/²⁰⁴Pb (x) and ²⁰⁷Pb/²⁰⁴Pb (y) ratios are needed to solve for the model age.

To apply CalcModAge() to a hypothetical sample with ²⁰⁶Pb/²⁰⁴Pb = 13.5 and ²⁰⁷Pb/²⁰⁴Pb = 14.5:

All calculations from here onwards require that the sample age is known. Preferably an independently obtained age, such as a zircon U–Pb age is used. When comparing Pb isotope signatures across different time periods, it is often more useful to compare the model source μ (²³⁸U/²⁰⁴Pb) rather than the Pb isotope ratios, as μ does not vary with time in a closed system.

This is implemented in PbIso by the function CalcMu(), using the sample age (t) in Ma, ²⁰⁶Pb/²⁰⁴Pb (x), and ²⁰⁷Pb/²⁰⁴Pb (y) ratios. Optional arguments T1, X1, Y1, U8U5, L5, and L8 can also be applied in the format CalcMu(t, x, y, T1, X1, Y1, U8U5, L5, L8). For example, if we have a sample with a known deposit age (independently obtained using a robust method such as U–Pb zircon crystallisation) of t = 2700 Ma; ²⁰⁶Pb/²⁰⁴Pb = 13.5, and ²⁰⁷Pb/²⁰⁴Pb = 14.5, and accepting the default Stacey and Kramers (1975) model values, we can calculate the model source μ by

These calculations can be visualised in Fig. 3b. The intersection of the isochron (m_sample; red line in Fig. 3b), the paleoisochron associated with the sample age (2700 Ma in this case), and the model source μ curve of 8.43 (blue curve in Fig. 3b) mark the initial Pb isotope composition of the sample. For samples with very low U concentrations, such as galena, the initial compositions will be approximately the same as the measured values. See Section 3.5 for calculating initial Pb isotope ratios, and Section 5.2 for generating the model curves and paleoisochron and isochron lines shown in Fig. 3b.

Somewhat similar to using the model source μ (²³⁸U/²⁰⁴Pb) for a sample, we can use the time-integrated κ (²³²Th/²³⁸U) to look at thorogenic Pb isotopic trends for samples or regions over different time scales.

This is implemented in PbIso as the CalcKa() function, using the sample age (t), ²⁰⁸Pb/²⁰⁴Pb (z) and ²⁰⁶Pb/²⁰⁴Pb (x). Again, let’s assume a sample with ²⁰⁶Pb/²⁰⁴Pb = 13.5; t = 2700 Ma; and now with ²⁰⁸Pb/²⁰⁴Pb = 33:

This calculation can be visualised in Fig. 4, where the red circle is the sample, the blue curve is evolution of the sample κ, and the grey circle is the initial Pb isotope ratio.

The above equations are identical to the Pb isotope equations in Section 3.1, so we could use those same functions, but substitute the sample μ (or ω) in. However, this would require two steps: (1) calculate the model source μ (or ω) using the CalcMu() (or CalcMu() ⋅ CalcKa() for ω), which requires t_sample, X_sample, Y_sample (and Z_sample for eq. 7) as input parameters and (2) use that calculated μ (or ω) value to input into equations Calc64(), Calc74(), and Calc84(). To eliminate the need to do this in two steps, we have added the functions Calc64in(), Calc74in(), and Calc84in() to calculate the initial Pb isotope ratios directly. Note that X_sample, Y_sample, and Z_sample are compulsory arguments because these are required to calculate μ and (or) κ. The initial Pb isotope ratios can be calculated using our hypothetical sample as follows:

Pb evolution models can be generated using the modelcurve() function in PbIso. This function takes the arguments for time (t)—usually given as a time interval e.g., 0:3700—as well as optional arguments for model starting parameters, μ₁, ω₁ values and decay constants (see Table 2). Additionally, parameters ϵ₁ and ϵ₂ can be specified to allow a variable μ with time (see Cumming and Richards (1975) for further information). These two ϵ parameters are rate factors that account for accelerated or decelerated Pb isotope evolution, and therefore the changes in μ of a Pb source over time. The modelcurve() function will generate a dataframe with columns t, x, y, and z. These correspond to time, ²⁰⁶Pb/²⁰⁴Pb, ²⁰⁷Pb/²⁰⁴Pb, and ²⁰⁸Pb/²⁰⁴Pb, respectively. The values x, y, and z are calculated following eqs. 2, 3, and 4 in Section 3.1, with the added parameters ϵ₁ and ϵ₂ to allow a variable μ through time. The model curves shown in Figs. 2c and 5 are generated using the modelcurve() function, and we have detailed the steps to do this below.

To generate a simple (Stacey and Kramers 1975) second stage Pb evolution model, only the time (t) is needed in the function modelcurve(). The “SKcurve” dataframe will have 3701 rows of data, each corresponding to a 1 Ma time interval (Table 4).

The modelcurve() function can be used to produce model curves with different model parameters. For example the “modelex1” curve in Fig. 5, which corresponds to a model μ of 8 is generated by the following:

Using the same parameters as above, but using a variable μ with time (specified by parameters ϵ₁ and ϵ₂)

The modelcurve() function can be used to produce model curves over different time ranges with optional arguments. For example, in the code below, “modelex3” is a hypothetical model using a custom model “my_model”. These three example curves are shown in Fig. 5.

Usually, in Pb isotope studies of ore deposits, multiple samples will be obtained from an individual deposit, which produce analyses with a range of Pb isotope values. When doing large regional-scale studies, often only the least radiogenic sample from each deposit will be used.

In PbIso, we have implemented the function LeastRad() that filters a dataset (e.g., df), based on the lowest analysis of an isotope ratio (e.g., ²⁰⁷Pb/²⁰⁴Pb or ²⁰⁶Pb/²⁰⁴Pb) from each group (e.g., ore deposit), in the format: LeastRad(df, group, value).

For example, to filter a sample dataset “df”, based on the lowest ²⁰⁷Pb/²⁰⁴Pb analysis from each deposit we would use:

Note that this function can only be applied to a dataframe in R, not to individual measurements. More information on using PbIso with dataframes is included in the sections below.

In addition to plotting sample data, there are also several plotting features such as paleoisochron and isochron lines. To generate paleoisochron lines for a given time (t), the slope and y-intercept are needed. To calculate the slope of a paleoisochron line on a ²⁰⁶Pb/²⁰⁴Pb versus ²⁰⁷Pb/²⁰⁴Pb plot, the function isochron76slope() is used, which takes the argument t as well as optional arguments (see documentation). The associated y-intercept for that paleoisochron is given by the function isochron76yint(). These can then be used to plot the paleoisochron line along with a model curve. Similarly, to calculate the paleoisochron slope and y-intercept on a ²⁰⁶Pb/²⁰⁴Pb versus ²⁰⁸Pb/²⁰⁴Pb plot, the functions isochron86slope() and isochron86yint() can be used. Paleoisochron lines (and isochron line for t = 0) in Fig. 2 are plotted using the below functions for t = 3000, 2000, 1000, and 0 Ma. These values can then be used to plot paleoisochron lines, or the function can be called directly, for example, by using abline(a = isochron76yint(2700), b = isochron76slope(2700)) in base R plotting, or geom_abline(slope = isochron76slope(2700), intercept = isochron76yint(2700)) in ggplot2 (this command is used to produce the black dashed line in 3.1b). For a ²⁰⁶Pb/²⁰⁴Pb versus ²⁰⁸Pb/²⁰⁴Pb plot, the ggplot2 command would be geom_abline(slope = isochron86slope(2700), intercept = isochron86yint(2700)) (black dashed line in 3.2). To use this function to generate the slope and y-intercept of a paleoisochron at time 2700 Ma:

We have shown above that the PbIso package allows for straightforward calculations of various Pb isotope parameters such as model age, model source μ (²³⁸U/²⁰⁴Pb), time-integrated κ (²³²Th/²³⁸U), and initial Pb isotope ratios. However, usually we will want to apply these calculations to an entire dataset rather than to just one sample. PbIso is packaged with a sample dataset, which is a subset of sulfide Pb isotope analyses from the Superior Province in Canada obtained from the DepIso database (see https://sil.usask.ca/databases.php and Eglington (2018a)). We briefly document below how to apply functions to this dataset, using the “SampleData.csv”.

Import the SampleData.csv file:

The dataframe that we have imported as “df”, before any calculations have been applied, is shown in Table 5.

We can now apply the PbIso functions to the dataframe in the same way we do to individual analyses. Each of the new calculations below will be added as separate columns to the “df” dataframe.

Like applying the functions to individual analyses, we can specify optional arguments (e.g., T1, X1, and Y1) or specify a predefined or user-defined model (see Section 2.2 for defining your own model “my_model”):

The model age function is slightly more complex, so we need to use the base R function mapply(). Instead of using “model” to add a predefined or user-defined model, we need to use “MoreArgs”.

The resulting “df” dataframe with calculations applied (we have removed the extra columns demonstrating the optional arguments and models) is shown in Table 6.

The dataframe is now ready to use for plotting various parameters against each other or for performing a wide range of statistical analyses that is possible with other R functions and packages.

Users may wish to only use the least radiogenic sample from each deposit, which can be performed using the LeastRad() function, either as the first step in this workflow (immediately after importing the dataset) or the last (after running the various calculations above).

In either case:

This produces a dataframe (dfLR) that contains only the sample with the lowest ²⁰⁷Pb/²⁰⁴Pb value from each deposit (DepositName).

Often Pb isotope models are developed to help understand the evolution of Pb isotopes in particular regions. Two models are commonly referred to when dealing with Superior Province data, the Abitibi–Wawa model (Thorpe 1999) and an Archean model based on sulfide data from the Pilbara Craton in Australia and other Archean terranes (Thorpe et al. 1992). We refer to these below as the “Abitibi Model” and “Archean Model”, respectively.

With PbIso, it is very straightforward to generate multistage models. In the hypothetical example below, let’s assume we want to model an extraction event from the Stacey and Kramers (1975) second stage model curve at 3200 Ma, with a new μ value of 5. First, the starting parameters need to be obtained. To do this, we can just filter the “SK2model” dataframe for our starting time, t = 3200 Ma, as follows:

The new_start_params values are now the starting composition for our “NewSuperior” model below, using our desired μ value of 5. We can then use the modelcurve() function to generate the dataframe for this model and plot it on Fig. 6.

The above steps can be repeated indefinitely to generate additional model “stages” using the PbIso functions, although caution should be applied to whether this is geologically reasonable. By plotting the Pb isotope data along with the model curves, we can begin to interrogate different Pb evolution models for what might be realistic for the source of Pb in sulphides from the Superior Province. Note that the “NewSuperior” model is very much a hypothetical example to demonstrate how this can be done in the PbIso package, and is not being suggested here as a suitable model for Pb isotope evolution in the Superior Province.

Using the above functions within R allows for flexibility in applying them to different datasets and using models that are appropriate for specific regions of interest. For a quick and easy to use approach, and to demonstrate some of the capabilities of PbIso, we have deployed the PbIso package into a Shiny application (see https://shereearmistead.github.io/software/pbiso). The app allows users to add their own data by simply copying and pasting data into an Excel-like spreadsheet (Fig. 7). The app will then automatically generate the values for model age, model source μ (²³⁸U/²⁰⁴Pb), time-integrated κ (²³²Th/²³⁸U), and the three initial Pb isotope ratios. The app also plots these data onto a series of standard plots. These include (1) model source μ (²³⁸U/²⁰⁴Pb) versus age; (2) time-integrated κ (²³²Th/²³⁸U) versus age; (3) ²⁰⁶Pb/²⁰⁴Pb versus ²⁰⁷Pb/²⁰⁴Pb; and (4) ²⁰⁶Pb/²⁰⁴Pb versus ²⁰⁸Pb/²⁰⁴Pb. The app also allows users to modify the model parameters such as T1, X1, Y1, and decay constants, or select from one of the included models.

The spreadsheet with newly calculated values can be exported as a .xlsx file by clicking the “Download data” button (Fig. 7). The downloaded spreadsheet contains two tabs, the “DataOutput” tab contains the input data with the new calculations as shown on the left-hand side of the Shiny App screen. The second tab in the .xlsx spreadsheet is the “ModelParameters” tab and contains all of the values from the “ModelParameters” tab in the Shiny App that were used to generate the calculated values. The four plots can be downloaded separately as .pdf files by clicking the “Download plot as pdf” button in the lower right corner (Fig. 8). If the input sample data are 12 samples or less, these will be differentiated by colour and the sample ID included in a legend below the plots. For more than 12 samples, differentiating by colour becomes difficult, and so these will simply be plotted as the same colour points (Fig. 8). Optional lines and curves for different Pb isotopic models can be selected or deselected using the controls in the lower panel. The x and y axis limits can be specified either by typing a number or using the up/down arrows in the y min, y max, x min, and x max fields. The two isotopic plots have the option of adding a 95% filled contour behind the data.

We have provided a user-friendly R package for dealing with Pb isotope data. The functions allow flexibility in that they can be used in a very simple way accepting the default values for the Stacey and Kramers (1975) second stage model, or the user can change individual parameters or apply a user-defined model. This toolset adds to the growing number of open-source software packages that help with processing and interpreting geological data. A preprint version of this manuscript is available through EarthArXiv (Armistead et al. 2021). This package may continue to have features added beyond the publication of this manuscript, and all updates will be managed and maintained through https://github.com/ShereeArmistead/PbIso.

Editor Fernando Corfu and reviewers Pieter Vermeesch, Tom Andersen, and an anonymous reviewer are thanked for their constructive feedback that improved the PbIso package and manuscript. Bradley Cave, Patrick Sack, Ryan Ickert, and Michael Anenburg are thanked for discussion and feedback that improved the R package and Shiny application.

All data used in the manuscript are included in the PbIso package https://github.com/ShereeArmistead/PbIso.

Name of the code/library: PbIso

Contact: [email protected]

Hardware requirements: …

Program language: R

Software required: R

Program size: 1 MB

The source codes are available for downloading at the link https://github.com/ShereeArmistead/PbIso.

Conceptualization: SEA, SJP

Data curation: SEA, BME

Formal analysis: SEA, BME

Funding acquisition: BME, SJP

Investigation: SEA, BME

Methodology: SEA, BME

Project administration: SEA, SJP

Resources: SJP

Software: SEA

Supervision: BME, SJP

Validation: SEA, SJP

Visualization: SEA

Writing – original draft: SEA, BME

Writing – review & editing: SEA, BME, SJP

SEA was supported by the Canada First Research Excellence Fund Metal Earth Programme at Laurentian University.