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

Abstract

ACCEPTED MANUSCRIPT Canadian Journal of Earth Sciences 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 206 Pb/ 204 Pb, 207 Pb/ 204 Pb or 208 Pb/ 204 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 238 U and 235 U produces very different decay characteristics for old (Archean) rocks and minerals when compared with younger samples.Since 235 U decays much more rapidly than 238 U, the 207 Pb/ 204 Pb ratio increases much more rapidly than 206 Pb/ 204 Pb in early Earth history.At younger (Mesoproterozoic to Phanerozoic) times, there is little change in 207 Pb/ 204 Pb while 206 Pb/ 204 Pb continues to increase (Figure 1).As with all isotope systems, complexity may be introduced by geological processes such as multiple dis-ACCEPTED MANUSCRIPT Canadian Journal of Earth Sciences turbance 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) which exclusively link to silicate petrogenetic processes, the U-Th-Pb isotope system provides insights for both silicate and sulphide processes.
Although not frequently used, the Re-Os isotope system also provides information on sulphide 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 'common Pb' isotope system was extensively used in the 1970's and 1980's, prior to development of ICPMS and ion microprobe capabilities for Lu-Hf investigations of zircon.It is often seen as a more difficult system to understand, in part because it uses graphs with radiogenic products on both axes, either 206 Pb/ 204 Pb vs 207 Pb/ 204 Pb for the uranogenic system or 206 Pb/ 204 Pb vs 208 Pb/ 204 Pb for the thorogenic system.As with other isotope decay systems, the isochron relationships for the U-Th-Pb system are described by relationships of the form: 206   204 Where $\lambda$ 8 is the decay constant for the radioactive isotope ( 238 U in this example), T initial is the age in years of the initial event and t is the age of a secondary event.For t=0 (present day) this latter term reduces to the value 1.The full set of equations for the three isotope systems are given in equations ( 2), ( 3) and (4).For simplicity, some ratios are often represented by Greek characters by convention. 238U/ 204 Pb is represented by the letter , 232 Th/ 204 Pb by the letter  and 232 Th/ 238 U by the letter .
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 and by Houtermans (Holmes (1946) and 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.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 (Eglington (2019), Johnson and DePaolo (1994)) or because the magmas involved tap material for each isotope system from multiple, distinct reservoirs (Stracke (2012), Vervoort et al. (1994)).
The viability of the latter has become increasingly important with recent suggestions that distinct zones of concentration of sulphide minerals may occur in the crust and mantle lithosphere which 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 (Be'eri-Shlevin et al. (2010), Carr et al. (1995), Champion and Huston (2016), Eglington (2018a), Gulson (1986), Halla (2018), Huston et al. (2016), Luais and Hawkesworth (2002), Warren and Thorpe (1994)).
Since 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 (Albarède et al. (2012), Andersen (1998), Champion and Huston (2016), Eglington (2018a), Gale and Mussett (1973), Halla (2018), Harmer et al. (1995)).The model source  ( 238 U/ 204 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), Eglington (2018b), Vermeesch (2018), Ludwig (2001)).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 (Faure (1977), Gale and Mussett (1973), Halla (2018)).

Introduction to PbIso
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 (GUI), 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 ACCEPTED MANUSCRIPT Canadian Journal of Earth Sciences 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.

Pb isotope systematics
The decay of  and Doe, 1981); Bulk Silicate Earth (Maltese and Mezger, 2020); or region specific models (e.g.Abitibi-Wawa in Superior Province (Thorpe, 1999)). 2 Using PbIso in R PbIso includes functions for calculating initial Pb isotope ratios, model age, model source  ( 238 U/ 204 Pb) and time-integrated  ( 232 Th/ 238 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

Installation
The package can be downloaded from https://github.com/ShereeArmistead/PbIsoor can be installed within R by running the following: install.packages(devtools) devtools::install_github("shereearmistead/PbIso") library(PbIso) 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.

Running functions
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:

Customising model parameters
As stated previously, the default starting parameters are based on Stacey and Kramers (1975) 2nd stage model, however, these can be manually overridden by specifying them in the function.For further information about how to generate a table with values for a customised Pb evolution model, see section 4. The optional parameters, in this case, T1, X1, Mu1, can be specified as: Calc64(2700, T1 = 4000, X1 = 9.5, Mu1 = 7) [1] 11.87765 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) 1st stage model: Calc64(2700, model = SK1) [1] 12.98544 Or the Maltese and Mezger (2020) Bulk Silicate Earth Model: Calc64(2700, model = MM20) [1] 13.57078 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:  Calc64(2700, model = my_model) [1] 13.89664 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.
Table 3: Predefined models that can be used in the PbIso functions.

Calc84(2700)
[1] 33.366 See Table 2 for the optional parameters for each of these functions to allow a different Pb evolution model to be used.

Model age
Pb-Pb model ages are calculated by assuming a starting composition (typically Stacey and Kramers (1975) 2nd stage values), and calculating the time needed to reach the present-day measured values.To visualise this, in Figure 3a, the line connecting the model starting composition (X1, Y1 at T1) and the sample (red circle) intersects the Stacey and Kramers (1975)  To apply CalcModAge() to a hypothetical sample with 206 Pb/ 204 Pb = 13.5 and 207 Pb/ 204 Pb = 14.5: CalcModAge(13.5, 14.5) [1] 2510  CalcMu(2700, 13.5, 14.5) [1] 8.43 These calculations can be visualised in Figure 3b.The intersection of the isochron (m sample ; red line in Figure 3b), the paleoisochron associated with the sample age (2700 Ma in this case), and the model source  curve of 8.43 (blue curve in Figure 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 Figure 3b.The time-integrated  ( 232 Th/ 238 U) is given by:

Time
This is implemented in PbIso as the CalcKa() CalcKa(2700, 33, 13.5) [1] 3.32 This calculation can be visualised in Figure 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.

Initial Pb isotope ratios
Often we are interested in the initial Pb isotope composition of a sample at the time of formation, particularly for more radiogenic samples.
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 equation ( 11)) 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: the 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: Calc64in(2700, 13.5, 14.5) [1] 13.304 Calc74in(2700, 13.5, 14.5) [1] 14.464 Calc84in(2700, 13.5, 14.5, 33) [1] 32.852

Pb evolution models
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,  1 ,  1 values and decay constants (see Table 2).Additionally, parameters  1 and  2 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, z.These correspond to time, 206 Pb/ 204 Pb, 207 Pb/ 204 Pb and 208 Pb/ 204 Pb, respectively.
The values x, y and z are calculated following equations ( 2), ( 3) and ( 4) in section 3.1, with the added parameters  1 and  2 to allow a variable  through time.The model curves shown in Figure 2c and Figure 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) 2nd 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).

Least radiogenic calculation
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() which filters a dataset (e.g.df), based on the lowest analysis of an isotope ratio (e.g. 207Pb/ 204 Pb or 206 Pb/ 204 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 207 Pb/ 204 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.

Plotting parameters
In addition to plotting sample data, there are also several plotting features such as paleoisochron and isochron lines.

Paleoisochrons
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 206 Pb/ 204 Pb vs. 207 Pb/ 204 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 isochron76slope(2700) [1] 0.681 isochron76yint(2700) [1] 5.4

Province, Canada
We have shown above that the PbIso package allows for straightforward calculations of various Pb isotope parameters such as model age, model source  ( 238 U/ 204 Pb), time-integrated  ( 232 Th/ 238 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.Using the standard base R function mapply() we can apply the PbIso functions to a dataframe.PbIso is packaged with a sample dataset, which is a subset of sulphide 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.csvfile: df <-read.csv("SampleData.csv") The dataframe that we have imported as 'df', before any calculations have been applied, is shown in Table 5.The model age function is slightly more complex, so we need to use the base R function mapply().Instead 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 207 Pb/ 204 Pb value from each deposit (DepositName).

Pb evolution models for Superior Province
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 sulphide 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.
AbitibiModel <-modelcurve(4000:0, model = THAW, Mu1 = 8) ArcheanModel <-modelcurve(4000:0, model = THAR, Mu1 = 9) With PbIso, it's very straightforward to generate multi-stage models.In the hypothetical example below, let's assume we want to model an extraction event from the Stacey and Kramers (1975)  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 Figure 6.
new_Superior_model <-list(T1 = 3200, X1 = 12.442, Y1 = 14.048,Z1 = 32.311,Mu1=5) new_Superior_curve <-modelcurve(3200:0, model = new_Superior_model) 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 geological 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.

Shiny application
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.is 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.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.

Conclusion
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) 2nd 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.

Figure 1 :
Figure 1: Single-stage growth curve.Evolution of the radiogenic isotopes 206 Pb and 207 Pb, relative to common 204 Pb, from t1 to the present-day.Reproduced from Halla (2018).

Figure 2 :
Figure 2: Pb isotopic evolution through time, a) evolution of 206 Pb/ 204 Pb and 207 Pb/ 204 Pb; b) evolution of 206 Pb/ 204 Pb and 208 Pb/ 204 Pb; both using the 2nd stage model parameters from Stacey and Kramers (1975).Paleoisochrons for 3000 Ma, 2000 Ma, 1000 Ma and the isochron for 0 Ma are shown as dashed lines in both plots.c) Pb evolution curves using models packaged within PbIso, see Table 3 for model parameters and references.Points are shown along each curve at 100 Ma intervals.See section 4 for how to generate the model curves and section 5.2.1 for the paleoisochron/isochron lines shown in (a) and (b).

Figure 3 :
Figure 3: Hypothetical sample (red circle) plotted on a standard Stacey and Kramers (1975) 2nd stage model curve (black line), showing a) that a line connected between a model starting composition and the sample composition, will project onto the Stacey and Kramers (1975) model curve at the corresponding model age, and b) additional values that can be calculated if the actual sample age (in this case 2700 Ma) is known.Grey circle represents the calculated initial ratios, the blue curve represents the corresponding model source  value for this sample, and the dashed line is the sample age (2700 Ma) paleoisochron.

Figure 4 :
Figure 4: Plot of 208 Pb/ 204 Pb vs 206 Pb/ 204 Pb with the same hypothetical sample as previous plots, with the 2700 Ma paleoisochron (dashed line), corresponding kappa curve (blue) and initial isotope composition (grey circle)

Figure 5 :
Figure 5: Model curves generated using the modelcurve() function for three hypothetical models.Points on curves are shown for every 100 Ma.

Figure 6
Figure 6: a) Stacey and Kramers (1975) model curve (black line) with three Pb evolution model curves, including two published models (Thorpe (1999) and Thorpe et al. (1992)), and a hypothetical NewSuperior model; and b) the same model curves as (a) but showing the extent of Superior Province Pb isotope data from selected ore deposits.Filled symbols are the least radiogenic value from each deposit and the unfilled symbols include all data from deposits.Circles along model curves are shown at 100 Ma intervals.

Figure 7 :
Figure 7: Screenshot of the PbIso shiny application.Left side is where the user can copy and paste data, and values are subsequently calculated.The right side has a series of tabs that include four commonly used Pb isotope plots and a tab for modifying the model parameters

Figure 8 :
Figure 8: PbIso Shiny app plots that are produced based on user input data 206 Pb/ 204 Pb or 207 Pb/ 204 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 238 U/ 204 Pb () which 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 238 U/ 204 Pb ratio needed 238U to 206 Pb; 235 U to 207 Pb and 232 Th to 208 Pb and the non-radiogenic 204 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 (Figure2), so samples with significantly different ages cannot easily be compared using the 206 Pb/ 204 Pb, 207 Pb/ 204 Pb, and 208 Pb/ 204Pb ratios.To allow easy comparison of samples across broad time periods, the model source  ( 238 U/ 204 Pb) and  ( 232 Th/ 238 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 Stacey and Kramers (1975)ser data including sample name/ID, age, 206 Pb/ 204 Pb, 207 Pb/ 204 Pb, and 208 Pb/ 204 Pb ratios.Users can then export the processed data as a .xlsxfile,whichwillinclude the calculated columns for initial Pb isotope ratios, model age, model source  ( 238 U/ 204 Pb) and time-integrated  ( 232 Th/ 238 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 .pdffigures.The functions in PbIso take one or more of the basic input parameters t (time (Ma)), x ( 206 Pb/ 204 Pb), y ( 207 Pb/ 204 Pb) and z ( 208 Pb/ 204 Pb) to perform the calculations.For advanced usage, the functions can also optionally take the values for different model parameters (summarised in Table1).The calculations and functions used in PbIso assume a starting composition and model followingStacey and Kramers (1975)2nd stage model, although this can easily be changed if an alternative model (e.g.Stacey and Kramers (1975) https://shereearmistead.github.io/software/pbiso),whichrequiresnoknowledge of R, making it accessible for users.singlestage;Malteseand Mezger (2020) Bulk Silicate Earth or others) is preferred.Table1: Default model parameters used in PbIso functions.These can be changed if an alternative model is preferred.

Table 2 :
Summary of 'PbIso' functions and their input parameters Paleoisochron slope on a 206 Pb/ 204 Pb vs 208 Pb/ 204 Pb plot

review of Pb isotopes and how to perform calculations with PbIso 3.1 The evolution of radiogenic Pb isotopes with time
Stacey and Kramers (1975)Pb) in equation (4) can also be expressed as  ⋅ , which is equivalent to 238 U/ 204 Pb ⋅ 232 Th/ 238 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.Forsimply calculating the value of each isotope ratio on theStacey and Kramers (1975)average ore lead curve at a given time, only the age is required.For example, the 206 Pb/ 204 Pb, 207 Pb/ 204 Pb and 208 Pb/ 204 Pb ratios on theStacey and Kramers (1975)curve at a given time, say 2700 Ma, is given by: curve at the model age for this sample.This is numerically given by the following equation: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 206 Pb/ 204 Pb (x) and 207 Pb/ 204 Pb (y) ratios are needed to solve for the model age.

3.3 Model source 𝜇 ( 238 U/ 204 Pb)
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  ( 238 U/ 204 Pb) rather than the Pb isotope ratios, as  does not vary with time in a closed system.Rocks or minerals that formed within the same reservoir will have Pb isotopic compositions that cluster along an isochron line on a 206 Pb/ 204 Pb vs 207 Pb 204 Pb plot (red line in Figure3b).The least radiogenic samples will fall near the lower left end of this isochron, while more radiogenic samples will fall near the upper right end. One wayto calculate a rock's age is to regress multiple sample analyses.The slope of this line is directly related to its mineralisation age.A more robust way of doing this though is to use the known age of a sample, whereby the slope (m) for a sample with a known age (t sample ) is defined by: Eglington (2018a)995)n PbIso by the function mslope(), which takes the argument t and additional optional arguments (see documentation).Substituting equation (6) into the following equation gives us the model source  ( 238 U/ 204 Pb)(Harmer et al. (1995),Eglington (2018a)): This is implemented in PbIso by the function CalcMu(), using the sample age (t) in Ma, 206 Pb/ 204 Pb (x)

-integrated 𝜅 ( 232 Th/ 238 U)
Somewhat similar to using the model source  ( 238 U/ 204 Pb) for a sample, we can use the time-integrated  ( 232 Th/ 238 U) to look at thorogenic Pb isotopic trends for samples or regions over different time scales.

Table 4 :
First ten rows of the 'SKcurve' dataframe produced using the modelcurve() function.

Table 5 :
First ten rows of sample data input with sample informa-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.

Table 6 :
First ten rows of dataframe after the PbIso functions have 2nd 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 app will then automatically generate the values for model age, model source  ( 238 U/ 204 Pb), time-integrated  ( 232 Th/ 238 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  ( 238 U/ 204 Pb) vs. age; 2) time-integrated  ( 232 Th/ 238 U) vs. age; 3) 206 Pb/ 204 Pb vs. 207 Pb/ 204 Pb; and 4) 206 Pb/ 204 Pb vs. 208 Pb/ 204 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.