Introduction to serac

Rosalie Bruel and Pierre Sabatier

Last update: 2024-11-15

The present vignette helps you compute your first age-depth model with serac, an R package for ShortlivEd RAdionuclide Chronology of recent sediment cores.

Note that the package evolves with time, so be sure to update your version when starting a new project.

Since its publication in 2020, we added a few options to the code. The most detailed list of arguments are in the help notice (?serac in the console). We are also including a few of the new options in this vignette (see sections [NEW] below)

1 First time user

1.1 Setting up the folders

Which folder are your currently working in?

getwd()

If you want to work in another folder, change it with setwd(), e.g.:

setwd("~/myProject/Age depth model")

Line below will create the Cores folder if it does not exist already

dir.create(file.path(getwd(), 'Cores'), showWarnings = FALSE)

Within the Cores folder, user needs to create one folder per core.

dir.create(file.path(paste0(getwd(),'/Cores'), 'serac_example_ALO09P12'), showWarnings = FALSE)

We are writing in this subfolder an example dataset (Source: <a href:“https://www.sciencedirect.com/science/article/pii/S0033589412000294”>Wilhelm et al., 2012).

Format your own data following this template, or get help with formatting using the serac_input_formatting() function.

?serac_example_ALO09P12
write.table(x = serac_example_ALO09P12, file = paste0(getwd(),'/Cores/serac_example_ALO09P12/serac_example_ALO09P12.txt'),col.names = T, row.names = F,sep="\t")
# Including proxy data for this core too
write.table(x = serac_example_ALO09P12_proxy, file = paste0(getwd(),'/Cores/serac_example_ALO09P12/serac_example_ALO09P12_proxy.txt'),col.names = T, row.names = F,sep="\t")

1.2 User metadata (optional)

One of our goal with serac is to improve reproducibility and data management.
A first and easy step is to associate the name of the scientist (you!) to each model built. To enter personal data (recommended: name, affiliation, ORCID number and email), tape in the following line in the console, and follow instructions.

user_infos()

2 Run serac

2.1 Minimum arguments

At the minimum, your function will be need the name of the core and the coring year. In the example below, I included the save_code argument because Rmarkdown does not save code history the way the code can use it. Unless you are including a plot in a Rmarkdown, I recommend you do not specify this argument, or turn it to its default, i.e., TRUE. That way, you can go back and see your code.

serac(name="serac_example_ALO09P12",coring_yr=2009,
      save_code=FALSE) #do not include 'save_code' if you are working outside a Rmarkdown document

2.2 Final model

Several sedimentation hypotheses were tested, and this is what the authors (Wilhelm et al, 2012) chose as the best model (type in only the function serac – everything that appears below (“Sediment accumulation rate …”) are model outputs):

model_ALO09P12 <-
  serac(name="serac_example_ALO09P12", coring_yr=2009, model=c("CFCS"),
      plotphoto=FALSE, minphoto=c(0), maxphoto=c(210),
      plot_Pb=T, plot_Am=T, plot_Cs=T, Cher=c(30,40), Hemisphere=c("NH"), NWT=c(51,61),
      sedchange=c(75.5),
      plot_Pb_inst_deposit=T, inst_deposit=c(20,28,100,107,135,142,158,186),
      suppdescriptor=TRUE, descriptor_lab=c("Ca/Fe"),
      historic_d=c(20,28,100,107,135,142,158,186),
      historic_a=c(1994,1920,1886,1868),
      historic_n=c("sept 1994 flood","1920 flood","1886 flood","1868 flood ?"), 
      min_yr=c(1750),
      dmax=c(180), 
      plotpdf=T, preview=T,
      save_code=FALSE)


 ________________________
Running serac for core 'serac_example_ALO09P12'.


 Sediment accumulation rate (CFCS model) 0-75.5mm: SAR= 1.18 mm/yr, R2= 0.93
                          Error:     +/- 0.1 mm/yr

 Sediment accumulation rate (CFCS model) 75.5mm-bottom: SAR= 0.8 mm/yr, R2= 0.923
                          Error:     +/- 0.08 mm/yr
 Approximation of age at change(s) in sedimentation rate:
     Best Age: 1952 (uncertainty: 1947-1957)


 Age approximation of instantaneous deposit(s) from CFCS model:
     The instantaneous deposit at 20-28 mm has an estimated range of: 1991-1994.
     The instantaneous deposit at 100-107 mm has an estimated range of: 1913-1929.
     The instantaneous deposit at 135-142 mm has an estimated range of: 1874-1898.
     The instantaneous deposit at 158-186 mm has an estimated range of: 1852-1880.

 Inventory (Lead): 10287.66 Bq/m2 (+/- 232.9 Bq/m2)
 Inventory (Cesium): 19023.688 Bq/m2 (+/- 149.7 Bq/m2)


 ________________________

 The calculation took 0.5153 secs.

If you want to see the preview, change preview=T in the code. Before that, make sure to extend your window in RStudio (large plot!)

2.3 Plot against mass_depth

The plotting options are quite flexible (colors, text size).
We also offer the possibility to plot the figure against mass depth, using mass accumulation rate instead of sediment accumulation rate for the CFCS model. The code is similar to what we used previously, you just need to add mass_depth = TRUE to your function.

model_ALO09P12_mass_depth <-
  serac(name="serac_example_ALO09P12", coring_yr=2009, preview = T,save_code=FALSE, Pbcol = "orange", modelcol = "midnightblue", mycex = .8,
      mass_depth=TRUE)


 ________________________
Running serac for core 'serac_example_ALO09P12'.


 Mass accumulation rate (CFCS model): MAR= 0.065 g/cm2/yr, R2= 0.959
                              Error:     +/- 0.003 g/cm2/yr
 Inventory (Lead): 11707.91 Bq/m2 (+/- 253.2 Bq/m2)
 Inventory (Cesium): 21549.328 Bq/m2 (+/- 151.9 Bq/m2)


 ________________________

 The calculation took 0.0244 secs.

Note that you can then decide to enter the depth data for the Cesium peaks, instantaneous deposit, and points to ignore in the unit g.cm-2. If you want to do that, add the argument input_depth_mm = FALSE.

2.4 Add varves counting

If varves counting are available, you can display the dates on the age-depth model.

To do so, include in your project’s folder a text file (tab delimited) with two columns: “Age” and “Depth”. Depths must be entered in mm. The file must be named with the name of your sediment core (e.g., “MyCore”), and the append “_varves.txt”.

For example, for the example core serac_example_ALO09P12, you would need to include in the serac_example_ALO09P12 folder a file named serac_example_ALO09P12_varves.txt. Below is an example of what the file would look like:

Age Depth
2009 0
2008 5
2007 10
2006 15
2005 20
2004 30
2003 35
2002 40
2001 45

2.5 Set detection limits for the radionuclides [NEW]

New arguments were added in March 2023 to enter detection limits for Pb, Cs and Am. The arguments are respectively DL_Pb, DL_Cs, and DL_Am. Any value below what is specified will be plotted with a different color and not included in the model (they are treated the same as the depths entered in the ignore argument).

An example of how to use it is provided below with the hiatus example.

If any value was below the detection limit, a message will be printed in the console

Example: 3 Pbex values were below the detection limit set by the user (25)

2.6 Add a hiatus [NEW]

As of November 15th, 2024, the function was modified to allow for hiatuses to be included in the model.

The description of the new argument and how to use it:

If there are any hiatus, enter it/them here. The input form is a list of as many vectors as there are hiatuses, and each vector is made of two elements minimum: the depth of the hiatus in mm and the length, in years. For example, enter ‘hiatus = list(c(100, 5))’ for a 5-years hiatus at 100mm. By default, “hiatus” will be printed on the graph. Optionally, the user can specify a name and a color to the by adding element to the vector in the list. Example: hiatus = list(c(100, 5, “hiatus”, “blue”), c(100, 5, ““,”red”)) will had two hiatuses. Only the first one will have a name, and they will each be plotted with a different colors. Default = NULL.

Below is an example of how to use the new argument hiatus, with the example data included in the R package. Note that the hiatuses are not real for this core and that mass depth is not the recommended approach.

model_ALO09P12_hiatus <-
  serac(name="serac_example_ALO09P12", coring_yr=2009, model=c("CFCS"),
        plotphoto=FALSE, minphoto=c(0), maxphoto=c(210), mass_depth = TRUE,
        plot_Pb=T, plot_Am=T, plot_Cs=T, Cher=c(30,40), Hemisphere=c("NH"), NWT=c(51,61),
        sedchange=c(75.5),
        plot_Pb_inst_deposit=T, inst_deposit=c(20,28,100,107,135,142,158,186),
        suppdescriptor=TRUE, descriptor_lab=c("Ca/Fe"),
        historic_d=c(20,28,100,107,135,142,158,186),
        historic_a=c(1994,1920,1886,1868),
        historic_n=c("sept 1994 flood","1920 flood","1886 flood","1868 flood ?"), 
        min_yr=c(1750),
        dmax=c(180), 
        plotpdf=T, preview=T,
        save_code=FALSE,
        plot_unit = "mm", 
        DL_Pb = 25, DL_Cs = 25, DL_Am = 3, 
        hiatus = list(c(50, 5, "test hiatus", "black"), c(132, 15)))


 ________________________
Running serac for core 'serac_example_ALO09P12'.

9 Cs values were below the detection limit set by the user (25): the Cs values and the corresponding Cs_er will be shown with a different colour on the plot.  
3 Pbex values were below the detection limit set by the user (25): the corresponding depths were added to the 'ignore' vector: 129.5, 140.5, 154.5.  

 Mass accumulation rate (CFCS model) 0-75.5mm: MAR= 0.065 g/cm2/yr, R2= 0.925
                              Error:     +/- 0.006 g/cm2/yr

 Mass accumulation rate (CFCS model) 75.5mm-bottom: MAR= 0.038 g/cm2/yr, R2= 0.895
                              Error:     +/- 0.005 g/cm2/yr
 Approximation of age at change(s) in sedimentation rate:
     Best Age: 1956 (uncertainty: 1951-1961)


 Age approximation of instantaneous deposit(s) from CFCS model:
     The instantaneous deposit at 20-28 mm has an estimated range of: 1997-1996.
     The instantaneous deposit at 100-107 mm has an estimated range of: 1919-1910.
     The instantaneous deposit at 135-142 mm has an estimated range of: 1861-1846.
     The instantaneous deposit at 158-186 mm has an estimated range of: 1838-1819.

 Inventory (Lead): 10227.167 Bq/m2 (+/- 231.2 Bq/m2)
 Inventory (Cesium): 19053.007 Bq/m2 (+/- 152.8 Bq/m2)


 ________________________

 The calculation took 0.0784 secs.

3 Explore outputs

Files with age-depth model and metadata are automatically saved in your working folder. If you assign your function to an object, serac will return a list with the parameters.

mymodel1 <- serac(name="serac_example_ALO09P12", coring_yr=2009, preview = T,save_code=FALSE)


 ________________________
Running serac for core 'serac_example_ALO09P12'.


 General message: It seems you already tried this code combination. 
 A historic of parameters tested can be looked up in the file
 'serac_model_history_serac_example_ALO09P12.txt' (in the core directory).
 ________________________
 


 Sediment accumulation rate (CFCS model): SAR= 1.2 mm/yr, R2= 0.96
                          Error:     +/- 0.05 mm/yr
 Inventory (Lead): 11707.91 Bq/m2 (+/- 253.2 Bq/m2)
 Inventory (Cesium): 21549.328 Bq/m2 (+/- 151.9 Bq/m2)


 ________________________

 The calculation took 0.0232 secs.

As always when computing a new code, serac will check whether you tried the code combination previously. Here, we used the same code than in the previous section, so the code display a message before the outputs.

Note that several objects are accessible from mymodel1:
- data,
- CFCS sediment accumulation rate,
- CFCS age-depth model,
- CFCS age-depth model interpolated,
- Inventories,
- plot.
Depending on the code you chose (and the number of models you visualize), the list of objects increases or decreases.

mymodel1$plot 

# Visualize the input data
knitr::kable(mymodel1$data, format = "markdown")
depth_min depth_max depth_top depth_bottom depth_avg mass_depth_top mass_depth_bottom mass_depth_avg depth_avg_event_corrected density Pb210ex Pb210ex_er Cs137 Cs137_er Am241 Am241_er Pbex Pbex_er Cs Cs_er Am Am_er thickness
0 5 0 5 2.5 0.000 0.235 0.1175 2.5 0.47 642.0 11.0 772.0 3.0 1.00 0.51 642.0 11.0 772.0 3.0 1.00 0.51 5
5 14 5 14 9.5 0.235 0.514 0.3745 9.5 0.31 504.0 14.0 535.0 4.3 1.30 0.90 504.0 14.0 535.0 4.3 1.30 0.90 9
14 21 14 21 17.5 0.514 0.780 0.6470 17.5 0.38 371.0 7.0 540.0 2.0 1.30 0.40 371.0 7.0 540.0 2.0 1.30 0.40 7
21 27 21 27 24.0 0.780 1.176 0.9780 24.0 0.66 217.0 6.0 446.0 2.0 1.60 0.40 217.0 6.0 446.0 2.0 1.60 0.40 6
27 34 27 34 30.5 1.176 1.680 1.4280 30.5 0.72 232.0 12.0 790.0 5.0 1.90 0.90 232.0 12.0 790.0 5.0 1.90 0.90 7
34 40 34 40 37.0 1.680 1.968 1.8240 37.0 0.48 305.0 11.0 1340.0 5.0 1.20 0.90 305.0 11.0 1340.0 5.0 1.20 0.90 6
40 45 40 45 42.5 1.968 2.258 2.1130 42.5 0.58 219.0 10.0 530.0 3.0 1.30 0.57 219.0 10.0 530.0 3.0 1.30 0.57 5
45 51 45 51 48.0 2.258 2.546 2.4020 48.0 0.48 206.0 10.0 416.0 3.3 3.20 0.70 206.0 10.0 416.0 3.3 3.20 0.70 6
51 56 51 56 53.5 2.546 2.846 2.6960 53.5 0.60 168.0 10.0 423.0 3.4 5.94 0.80 168.0 10.0 423.0 3.4 5.94 0.80 5
56 61 56 61 58.5 2.846 3.136 2.9910 58.5 0.58 147.0 8.0 395.0 4.0 5.40 0.70 147.0 8.0 395.0 4.0 5.40 0.70 5
61 67 61 67 64.0 3.136 3.406 3.2710 64.0 0.45 116.0 5.0 209.0 1.6 2.30 0.40 116.0 5.0 209.0 1.6 2.30 0.40 6
67 72 67 72 69.5 3.406 3.701 3.5535 69.5 0.59 122.0 9.0 181.0 2.3 2.60 0.70 122.0 9.0 181.0 2.3 2.60 0.70 5
72 79 72 79 75.5 3.701 4.016 3.8585 75.5 0.45 125.0 8.0 117.0 1.8 1.50 0.70 125.0 8.0 117.0 1.8 1.50 0.70 7
79 83 79 83 81.0 4.016 4.268 4.1420 81.0 0.63 113.0 8.0 62.4 1.3 1.35 0.60 113.0 8.0 62.4 1.3 1.35 0.60 4
83 89 83 89 86.0 4.268 4.520 4.3940 86.0 0.42 122.0 5.0 34.7 0.6 0.71 0.40 122.0 5.0 34.7 0.6 0.71 0.40 6
89 96 89 96 92.5 4.520 4.856 4.6880 92.5 0.48 58.8 3.0 16.2 0.3 NA NA 58.8 3.0 16.2 0.3 NA NA 7
96 102 96 102 99.0 4.856 5.198 5.0270 99.0 0.57 48.4 2.0 4.5 0.1 NA NA 48.4 2.0 4.5 0.1 NA NA 6
102 107 102 107 104.5 5.198 5.638 5.4180 104.5 0.88 29.6 4.0 5.2 0.3 NA NA 29.6 4.0 5.2 0.3 NA NA 5
107 113 107 113 110.0 5.638 6.010 5.8240 110.0 0.62 26.8 3.0 6.1 0.2 NA NA 26.8 3.0 6.1 0.2 NA NA 6
113 121 113 121 117.0 6.010 6.338 6.1740 117.0 0.41 27.4 3.0 9.4 0.2 NA NA 27.4 3.0 9.4 0.2 NA NA 8
121 127 121 127 124.0 6.338 6.632 6.4850 124.0 0.49 29.0 4.1 11.3 0.5 NA NA 29.0 4.1 11.3 0.5 NA NA 6
127 132 127 132 129.5 6.632 6.957 6.7945 129.5 0.65 23.0 3.2 12.1 0.3 NA NA 23.0 3.2 12.1 0.3 NA NA 5
132 149 132 149 140.5 6.957 8.011 7.4840 140.5 0.62 15.4 3.0 7.7 0.3 NA NA 15.4 3.0 7.7 0.3 NA NA 17
149 160 149 160 154.5 8.011 8.704 8.3575 154.5 0.63 11.0 2.0 5.5 0.2 NA NA 11.0 2.0 5.5 0.2 NA NA 11
# Visualize the inventories (if density was in the input data file)
mymodel1$Inventories
              Inventory      min      max
Lead   Bq.m-2  11707.91 11454.75 11961.07
Cesium Bq.m-2  21549.33 21397.45 21701.20
# Save the plot with custom parameters, e.g., Courier font and grey scale without editing default colors.
pdf(paste0(getwd(),"/Cores/serac_example_ALO09P12/serac_example_ALO09P12_custom.pdf"),width = 10, height = 5, family="Courier", colormodel = "grey")
mymodel1$plot
dev.off()
quartz_off_screen 
                2 

4 Other functions

4.1 help_serac()

Another way to read some help.

4.2 serac_input_formatting()

If you unsure about the input format of your core (if you have depth in mm, or column names that do not follow our example), you can use the function serac_input_formatting(name = “MyCore”) (MyCore being the name of your core and folder in which the data are located) and follow the steps.

4.3 core_metadata()

Enter here a more extensive list of parameters associated to your core. Current list of parameters includes:

4.4 serac_map()

If you entered GPS coordinates for your core in the previous step, tape in serac_map() into your console and view the location of all the systems you entered GPS coordinates for.
If you enter a vector, e.g., serac_map(c(“Lake1”, “Lake2”, “Lagoon1”)), then the map will be generated solely for these systems.

5 Citation

If you use serac, please cite our paper describing the package:

Bruel, R., Sabatier, P., 2020. serac: an R package for ShortlivEd RAdionuclide chronology of recent sediment cores. Journal of Environmental Radioactivity 225, 106449. https://doi.org/10.1016/j.jenvrad.2020.106449