The present vignette helps you compute your first age-depth model with serac, a R package for ShortlivEd RAdionuclide Chronology of recent sediment cores.
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 doesn’t 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_ALO09P12write.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")
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()
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 doesn’t save code history the way the code can use it. Unless you’re 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're working outside a Rmarkdown document
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)
rate (CFCS model) 0-75.5mm: SAR= 1.18 mm/yr, R2= 0.93
Sediment accumulation : +/- 0.1 mm/yr
Error
rate (CFCS model) 75.5mm-bottom: SAR= 0.8 mm/yr, R2= 0.923
Sediment accumulation : +/- 0.08 mm/yr
Errorchange(s) in sedimentation rate:
Approximation of age at : 1952 (uncertainty: 1947-1957)
Best Age
deposit(s) from CFCS model:
Age approximation of instantaneous 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.
The instantaneous deposit at
Inventory (Lead): 10287.66 Bq/m2 (+/- 232.9 Bq/m2)
Inventory (Cesium): 19023.688 Bq/m2 (+/- 149.7 Bq/m2)
________________________
2.5029 secs. The calculation took
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!)
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)
rate (CFCS model): MAR= 0.065 g/cm2/yr, R2= 0.959
Mass accumulation : +/- 0.003 g/cm2/yr
ErrorInventory (Lead): 11707.91 Bq/m2 (+/- 253.2 Bq/m2)
Inventory (Cesium): 21549.328 Bq/m2 (+/- 151.9 Bq/m2)
________________________
0.0884 secs. The calculation took
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.
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 |
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.
<- serac(name="serac_example_ALO09P12", coring_yr=2009, preview = T,save_code=FALSE)
mymodel1
: It seems you already tried this code combination.
General messagein the file
A historic of parameters tested can be looked up 'serac_model_history_serac_example_ALO09P12.txt' (in the core directory).
________________________
rate (CFCS model): SAR= 1.2 mm/yr, R2= 0.96
Sediment accumulation : +/- 0.05 mm/yr
ErrorInventory (Lead): 11707.91 Bq/m2 (+/- 253.2 Bq/m2)
Inventory (Cesium): 21549.328 Bq/m2 (+/- 151.9 Bq/m2)
________________________
0.0625 secs. The calculation took
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.
$plot mymodel1
# Visualize the input data
::kable(mymodel1$data, format = "markdown") knitr
depth_min | depth_max | density | Pb210ex | Pb210ex_er | Cs137 | Cs137_er | Am241 | Am241_er | Pbex | Pbex_er | Cs | Cs_er | Am | Am_er | depth_top | depth_bottom | depth_avg | thickness | depth_avg_2 | d |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 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 | 0 | 5 | 2.5 | 5 | 2.5 | 2.5 |
5 | 14 | 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 | 5 | 14 | 9.5 | 9 | 9.5 | 9.5 |
14 | 21 | 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 | 14 | 21 | 17.5 | 7 | 17.5 | 17.5 |
21 | 27 | 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 | 21 | 27 | 24.0 | 6 | 24.0 | 24.0 |
27 | 34 | 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 | 27 | 34 | 30.5 | 7 | 30.5 | 30.5 |
34 | 40 | 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 | 34 | 40 | 37.0 | 6 | 37.0 | 37.0 |
40 | 45 | 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 | 40 | 45 | 42.5 | 5 | 42.5 | 42.5 |
45 | 51 | 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 | 45 | 51 | 48.0 | 6 | 48.0 | 48.0 |
51 | 56 | 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 | 51 | 56 | 53.5 | 5 | 53.5 | 53.5 |
56 | 61 | 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 | 56 | 61 | 58.5 | 5 | 58.5 | 58.5 |
61 | 67 | 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 | 61 | 67 | 64.0 | 6 | 64.0 | 64.0 |
67 | 72 | 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 | 67 | 72 | 69.5 | 5 | 69.5 | 69.5 |
72 | 79 | 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 | 72 | 79 | 75.5 | 7 | 75.5 | 75.5 |
79 | 83 | 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 | 79 | 83 | 81.0 | 4 | 81.0 | 81.0 |
83 | 89 | 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 | 83 | 89 | 86.0 | 6 | 86.0 | 86.0 |
89 | 96 | 0.48 | 58.8 | 3.0 | 16.2 | 0.3 | NA | NA | 58.8 | 3.0 | 16.2 | 0.3 | NA | NA | 89 | 96 | 92.5 | 7 | 92.5 | 92.5 |
96 | 102 | 0.57 | 48.4 | 2.0 | 4.5 | 0.1 | NA | NA | 48.4 | 2.0 | 4.5 | 0.1 | NA | NA | 96 | 102 | 99.0 | 6 | 99.0 | 99.0 |
102 | 107 | 0.88 | 29.6 | 4.0 | 5.2 | 0.3 | NA | NA | 29.6 | 4.0 | 5.2 | 0.3 | NA | NA | 102 | 107 | 104.5 | 5 | 104.5 | 104.5 |
107 | 113 | 0.62 | 26.8 | 3.0 | 6.1 | 0.2 | NA | NA | 26.8 | 3.0 | 6.1 | 0.2 | NA | NA | 107 | 113 | 110.0 | 6 | 110.0 | 110.0 |
113 | 121 | 0.41 | 27.4 | 3.0 | 9.4 | 0.2 | NA | NA | 27.4 | 3.0 | 9.4 | 0.2 | NA | NA | 113 | 121 | 117.0 | 8 | 117.0 | 117.0 |
121 | 127 | 0.49 | 29.0 | 4.1 | 11.3 | 0.5 | NA | NA | 29.0 | 4.1 | 11.3 | 0.5 | NA | NA | 121 | 127 | 124.0 | 6 | 124.0 | 124.0 |
127 | 132 | 0.65 | 23.0 | 3.2 | 12.1 | 0.3 | NA | NA | 23.0 | 3.2 | 12.1 | 0.3 | NA | NA | 127 | 132 | 129.5 | 5 | 129.5 | 129.5 |
132 | 149 | 0.62 | 15.4 | 3.0 | 7.7 | 0.3 | NA | NA | 15.4 | 3.0 | 7.7 | 0.3 | NA | NA | 132 | 149 | 140.5 | 17 | 140.5 | 140.5 |
149 | 160 | 0.63 | 11.0 | 2.0 | 5.5 | 0.2 | NA | NA | 11.0 | 2.0 | 5.5 | 0.2 | NA | NA | 149 | 160 | 154.5 | 11 | 154.5 | 154.5 |
# Visualize the inventories (if density was in the input data file)
$Inventories
mymodel1
Inventory min max-2 11707.91 11454.75 11961.07
Lead Bq.m-2 21549.33 21397.45 21701.20
Cesium Bq.m# 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")
$plot
mymodel1dev.off()
png 2
Another way to read some help.
If you unsure about the input format of your core (if you have depth in mm, or column names that don’t 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.
Enter here a more extensive list of parameters associated to your
core. Current list of parameters includes:
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.
If you use serac, please cite our paper describing the package:
Bruel, R., Sabatier, P., 2020. serac: A 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