R/SuessR_Package.R
In ctclark/SuessR: Suess and Laws Corrections for Marine Stable Carbon Isotope Data
Defines functions
reg.uptake
SuessR.custom
SuessR
laws.fun
Documented in
reg.uptake
SuessR
SuessR.custom
# SuessR Package
globalVariables("SuessR.reference.data")
# Function to calculate the Laws correction. Used within the SuessR() and SuessR.custom() functions.
laws.fun <- function(e1, e2, e.1, laws.CO2, sst, r, b, p) {
e2 + e1 - e.1 - (1/(1+laws.CO2*((((p)*(86400)*4*pi*(r*10^-6)^2))/((0.89*1.47^((sst + 273.15 - 303.15)/10))*((0.216*((4/3)*pi*r^3)^0.939)*10^-12)*(1+b)))))*((e2 - e.1)/(b+1))
}
# Function to calculate and apply Suess and Laws corrections.
#' @encoding UTF-8
#' @title Calculate Suess and Laws corrections for stable carbon isotope data from a built-in region
#' @description Generates region-specific Suess, Laws, and net (Suess + Laws) corrections for stable carbon isotope data
#' \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data input by the user. The net correction is then used to calculate the corrected \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data, which are supplied
#' in the output. This function is specifically for data from regions currently built into SuessR
#' ("Bering", "Aleutians","Gulf of Alaska", and "Subpolar North Atlantic", as of Version 0.1.3).
#' \if{html}{\figure{Built-in-regions.png}{options: width=600 alt="Built-in regions"}}
#' \if{latex}{\figure{Built-in-regions.png}{options: width=6in}}
#' @param data A matrix or data frame including columns containing sample ID ('id'), year of sample collection ('year'),
#' region of sample collection ('region'), and uncorrected \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C
#' data ('d13c').
#' @param correct.to The year to which the \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data will be corrected. Defaults to 1850 (Industrial Revolution).
#' @details The SuessR() function uses the year and region of sample collection to calculate Suess, Laws, and net (Suess + Laws)
#' corrections for \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data from marine organisms. The Suess Correction represents the change in \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C values of dissolved
#' inorganic carbon (DIC) in the surface ocean, and is calculated using an exponential function, calibrated to the global decline
#' DIC \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C values. Surface waters in different regions exhibit varying rates of CO\ifelse{html}{\out{<sub>2</sub>}}{\eqn{_2}} uptake from the atmosphere as a result of
#' differing water mass properties and residence time at the surface. Thus, a region-specific modifier has been
#' calculated from empirical observations of changes in DIC \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C values through time and is applied to each region. The Laws
#' correction accounts for changes in stable carbon isotope fractionation during CO\ifelse{html}{\out{<sub>2</sub>}}{\eqn{_2}} uptake and photosynthesis by phytoplankton,
#' and is impacted by changes in aqueous CO\ifelse{html}{\out{<sub>2</sub>}}{\eqn{_2}} concentrations, temperature, and salinity, as well as community growth rates,
#' average cell diameter, average organic carbon content of phytoplankton cells, permeability of the phytoplankton plasmalemma to CO\ifelse{html}{\out{<sub>2</sub>}}{\eqn{_2}},
#' and the ratio of net diffusional loss of CO\ifelse{html}{\out{<sub>2</sub>}}{\eqn{_2}} to carbon fixation. Historic observations and reconstructions of atmospheric CO\ifelse{html}{\out{<sub>2</sub>}}{\eqn{_2}} concentrations,
#' sea surface temperature, and sea surface salinity are used to make these calculations. See references for more details.
#'
#' This function allows users to correct their data to a specific year, using the 'correct.to' argument. This argument defaults
#' to AD1850, which represents onset of the Industrial Revolution and the initiation of the large-scale combustion of fossil
#' fuels that created the Suess Effect. Correcting data to 1850 allows them to be compared to data from any samples collected during
#' or before 1850 (including archaeological samples), as well as to any other samples Suess corrected to the year 1850. Users
#' examining contemporary datasets may wish to correct their data to another year. For example, a user examining a dataset
#' of \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data spanning the years 1970-2010 might choose to correct all their data back to the year 1970, to correct the
#' older samples forward to 2010, or to correct all the samples to 1990, the middle of the time series. In any of these scenarios,
#' corrected data would be comparable to one another. When choosing a value for 'correct.to' consideration should be given to
#' the comparability of the user's data to existing data sets, clarity of presentation of results (i.e., clear statements of the
#' year to which the data were corrected), and reproducibility of results (i.e., presentation of both uncorrected and corrected
#' data so future users can repeat analyses or correct the raw data to another year).
#' @return The output of this function is a data frame that includes the sample ID ('id'), year ('year'),
#' uncorrected \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C ('d13c.uncor'), Laws correction ('Laws.cor'), Suess Correction ('Suess.cor'), net correction
#' ('net.cor'), and corrected \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C ('d13c.cor') for each sample. The corrected \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data is equal to the uncorrected \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data
#' plus the net correction. The units for all values are the standard 'per mil' used for \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data.
#' @examples
#' example.data <- data.frame(id = c("Sample 1", "Sample 2", "Sample 3", "Sample 4",
#' "Sample 5", "Sample 6", "Sample 7", "Sample 8"),
#' year = c(2017, 2017, 2017, 2017, 1977, 1977, 1977, 1977),
#' d13c = c(-12, -12, -12, -12, -12, -12, -12, -12),
#' region = c("Bering Sea", "Aleutian Islands",
#' "Gulf of Alaska", "Subpolar North Atlantic",
#' "Bering Sea", "Aleutian Islands",
#' "Gulf of Alaska", "Subpolar North Atlantic"))
#'
#' SuessR(data = example.data)
#' @references Clark, C.T., M.R. Cape, M.D. Shapley, F.J. Mueter, B.P. Finney, and N. Misarti. (In Prep) SuessR: Regional Suess
#' and Laws corrections for \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C from marine organisms.
#'
#' Clark, C.T., L. Horstmann, A. de Vernal, A.M. Jensen, and N. Misarti. (2019) Pacific walrus diet across 4000 years of
#' changing sea ice conditions. \emph{Quaternary Research}, 1-17.
#'
#' Misarti, N., B. Finney, H. Maschner, and M.J. Wooller. (2009) Changes in northeast Pacific marine ecosystems over the last 4500
#' years: evidence from stable isotope analysis of bone collagen from archaeological middens. \emph{The Holocene}, 19:8. 1139-1151.
#'
#' @importFrom stats na.omit
#' @export
SuessR <- function(data, correct.to = 1850) {
# Make sure the data is a dataframe
if(is.matrix(data)) data <- as.data.frame(data)
# Make input data column names are correct
if(any(!(colnames(data) %in% c("id", "year", "region", "d13c")))) stop(paste("Unrecognized input data:", unique(colnames(data)[!(colnames(data) %in% c("id", "year", "region", "d13c"))]), " "))
# Screen input data to check for (and remove) missing values, with warning
if(anyNA(data)) warning("Rows with missing values were deleted")
data <- na.omit(data)
# Screen input data to make sure all years provided are integers between 1850 and the current year -1
current.year <- as.numeric(substr(date(), 21,24))-1
if(any(!(data$year %in% 1850:current.year))) stop(paste("Year must be an integer between 1850 and ", current.year))
# Screen input data to make sure the 'correct.to' year is an integer between 1850 and the current year -1
current.year <- as.numeric(substr(date(), 21,24))-1
if(any(!(correct.to %in% 1850:current.year))) stop(paste("correct.to must be an integer between 1850 and ", current.year))
# Screen input data to check for non-negative d13c values
if(any(data$d13c >= 0)) warning(paste("Some d13c value(s) are non-negative"))
# Screen input data to check for (and remove) missing values, with warning
if(any(!(data$region %in% SuessR.reference.data$region))) stop(paste("Unrecognized region:", unique(data$region[!(data$region %in% SuessR.reference.data$region)]), " "))
# Create output data frame
SuessR.out <- data.frame(id = data$id, year = data$year)
ref <- SuessR.reference.data
#ln(K0)
ref$lnK0 <- with(ref,
-58.0931 + 90.5069 * (100/(sst + 273.15))
+ (22.294 * log((sst + 273.15)/100))
+ (sss * (0.027766 + -0.025888 * ((sst+273.15)/100)
+ 0.005058 * ((sst + 273.15)/100)^2)))
#Ocean Increase
ref$oi <- NA
for(i in unique(ref$region)) {
sub.i <- ref[ref$region == i,]
n <- nrow(sub.i)
ref[ref$region == i, ]$oi <- c(NA, with(sub.i, (CO2atm[-1] - CO2atm[-n]) * sub.i$Cp[1]))
}
#~f(CO2)ocean
ref$fCO2 <- NA
for (i in unique(ref$region)) {
ref[ref$region == i,]$fCO2 <- c(285.78, cumsum(ref[ref$region == i, ]$oi[-1]) + 285.78)
}
#CO2aq
ref$CO2aq <- exp(ref$lnK0)*ref$fCO2
# Laws expression for a given year
ref$laws.current <- with(ref, laws.fun(e1=1, e2=26.5, e.1=1, p=1.5e-5, laws.CO2=CO2aq, sst=sst, r=r, b=0.2))
data$order <- seq(1,length(data$id),1)
data <- merge(data, ref, c("region", "year"), sort = F)
dat1850 <- ref[ref$year==1850,]
dat1850$laws1850 <- with(dat1850, laws.fun(e1=1, e2=26.5, e.1=1, p=1.5e-5, laws.CO2=CO2aq, sst=sst, r=r, b=0.2))
data <- merge(data, dat1850[,c("region", "laws1850")], "region", sort = F)
dat.correct.to <- ref[ref$year==correct.to,]
dat.correct.to$laws.correct.to <- with(dat.correct.to, laws.fun(e1=1, e2=26.5, e.1=1, p=1.5e-5, laws.CO2=CO2aq, sst=sst, r=r, b=0.2))
data <- merge(data, dat.correct.to[,c("region", "laws.correct.to")], "region", sort = F)
data$Laws.cor <- with(data, (laws.current - laws1850) - (laws.correct.to - laws1850))
SuessR.out <- data[,c("id","year", "d13c", "Laws.cor")]
names(SuessR.out)[3] <- "d13c.uncor"
SuessR.out$Suess.cor <- with(data, up.con*exp((year-1850)*0.027)
- up.con*exp((correct.to-1850)*0.027))
SuessR.out$net.cor <- SuessR.out$Suess.cor + SuessR.out$Laws.cor
SuessR.out$d13c.cor <- data$d13c + SuessR.out$net.cor
SuessR.out <- SuessR.out[order(data$order),]
print(cbind(SuessR.out[,1:3], round(SuessR.out[,4:7],3)))
}
# A function that allows users to add custom regions for the SuessR calculations
#' @encoding UTF-8
#' @title Calculate Suess and Laws corrections for stable carbon isotope data from a custom region
#' @description The SuessR.custom() function generates region-specific Suess, Laws, and net (Suess + Laws) corrections for
#' \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data input by the user. The net correction is then used to calculate the corrected \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data,
#' which are supplied in the output. This function is specifically for data from regions not currently built into SuessR
#' (i.e., different from "Bering", "Aleutians", "Gulf of Alaska", and "Subpolar North Atlantic", as of Version 0.1.3). Or
#' for scenarios in which users want to change the default parameters used by SuessR for the built-in regions (e.g., averge phytoplankton cell radius).
#' @param data A dataframe including sample ID, year of sample collection, uncorrected \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C
#' data, and region. This function is specifically for data from regions not
#' currently built into SuessR. Columns must be named 'id', 'year', 'd13c', and 'region'.
#' @param custom.region.data A data frame containing environmental data for the custom region from which the
#' samples originated. Must contain columns titled 'year', 'region', 'r', 'sst', 'sss', 'CO2atm', 'up.con', and 'Cp'.
#' See details for information on how to supply these parameters appropriately.
#' @param correct.to The year to which the \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data will be corrected. Defaults to 1850, pre-Suess effect.
#' @details The SuessR.custom() allows users to calculate and apply Suess and Laws corrections to \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data from marine
#' organisms collected in a region not currently built into the SuessR package. In Version 0.1.3, the
#' built-in regions are the Bering Sea ('Bering'), the Aleutian archipelago ('Aleutians'), and the Gulf of Alaska ('Gulf of
#' Alaska'). Because the Suess and Laws corrections require region-specific environmental data from 1850-present, users must
#' supply these data using the 'custom.region.data' argument. Once these data have been supplied, this function calculates
#' the corrections exactly like the SuessR() function. See the built-in 'SuessR.reference.data' object for an example
#' template for the data, as well as the references for detailed information on how to compile and supply the appropriate
#' data to the SuessR.custom() function. After compiling these data, please consider sharing them with the package authors
#' (ctclark 'at' alaska.edu) to be included as built-in regions in future version of this package.
#'
#' As with the SuessR() function, SuessR.custom() allows users to correct their data to a specific year, using the 'correct.to'
#' argument. This argument defaults to AD1850, which represents onset of the Industrial Revolution and the initiation of the
#' large-scale combusion of fossil fuels that created the Suess Effect. Correcting data to 1850 allows them to be compared to
#' data from any samples collected during or before 1850 (including archaeological samples), as well as to any other samples
#' Suess corrected to the year 1850. Users examining contemporary datasets may wish to correct their data to another time year.
#' For example, a user examining a dataset of \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data spanning the years 1970-2010 might choose to correct all their data back
#' to the year 1970, to correct the older samples forward to 2010, or to correct all the samples to 1990, the middle of the time
#' series. In any of these scenarios, corrected data would be comparable to one another. When choosing a value for 'correct.to'
#' consideration should be given to the comparability of the user's data to existing data sets, clarity of presentation of
#' results (i.e., clear statements of the year to which the data were corrected), and reproducibility of results (i.e.,
#' presentation of both uncorrected and corrected data so future users can repeat analyses or correct the raw data to another year).
#' @return The output of this function is a data frame that includes the sample ID ('id'), year ('year'),
#' uncorrected \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C ('d13c.uncor'), Laws correction ('Laws.cor'), Suess Correction ('Suess.cor'), net correction
#' ('net.cor'), and corrected \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C ('d13c.cor') for each sample. The corrected \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data is equal to the uncorrected \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data
#' plus the net correction. The units for all values are the standard 'per mil' used for \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data.
#' @examples
#' example.region.data <- data.frame(year = seq(from = 1850, to = 2020, by = 1),
#' region = rep("Example Region", 171),
#' r = rep(5, 171),
#' sst = seq(5.9, 6.6, 171),
#' sss = seq(32.3, 32.7, 171),
#' CO2atm = SuessR.reference.data$CO2atm,
#' up.con = rep(0.014, 171),
#' Cp = rep(0.41, 171))
#'
#' example.custom.data <- data.frame(id = c("Sample 1", "Sample 2",
#' "Sample 3", "Sample 4"),
#' year = c(1850, 1900, 1950, 2000),
#' d13c = c(-12, -12, -12, -12),
#' region = rep("Example Region", 4))
#'
#' SuessR.custom(data = example.custom.data,
#' custom.region.data = example.region.data,
#' correct.to = 1850)
#'
#' @importFrom stats na.omit
#' @export
SuessR.custom <- function(data, custom.region.data, correct.to = 1850) {
# Make sure the data is a dataframe
if(is.matrix(data)) data <- as.data.frame(data)
# Screen input data to check for (and remove) missing values, with warning
if(anyNA(data)) warning("Rows with missing values were deleted")
data <- na.omit(data)
# Make input data column names are correct
if(any(!(colnames(data) %in% c("id", "year", "region", "d13c")))) stop(paste("Unrecognized input data:", unique(colnames(data)[!(colnames(data) %in% c("id", "year", "region", "d13c"))]), " "))
# Screen input data to make sure all years provided are integers between 1850 and the current year -1
current.year <- as.numeric(substr(date(), 21,24))-1
if(any(!(data$year %in% 1850:current.year))) stop(paste("Year must be an integer between 1850 and ", current.year))
# Screen input data to make sure the 'correct.to' year is an integer between 1850 and the current year -1
current.year <- as.numeric(substr(date(), 21,24))-1
if(any(!(correct.to %in% 1850:current.year))) stop(paste("correct.to must be an integer between 1850 and ", current.year))
# Screen input data to check for non-negative d13c values
if(any(data$d13c >= 0)) warning(paste("Some d13c value(s) are non-negative"))
# Create output data frame
SuessR.out <- data.frame(id = data$id, year = data$year)
ref <- rbind(SuessR.reference.data, custom.region.data)
# Screen input data to check for (and remove) missing values, with warning
if(any(!(data$region %in% ref$region))) stop(paste("Unrecognized region: "), data$region[!(data$region %in% ref$region)])
#ln(K0)
ref$lnK0 <- with(ref,
-58.0931 + 90.5069 * (100/(sst + 273.15))
+ (22.294 * log((sst + 273.15)/100))
+ (sss * (0.027766 + -0.025888 * ((sst+273.15)/100)
+ 0.005058 * ((sst + 273.15)/100)^2)))
#Ocean Increase
ref$oi <- NA
for(i in unique(ref$region)) {
sub.i <- ref[ref$region == i,]
n <- nrow(sub.i)
ref[ref$region == i, ]$oi <- c(NA, with(sub.i, (CO2atm[-1] - CO2atm[-n]) * sub.i$Cp[1]))
}
#~f(CO2)ocean
ref$fCO2 <- NA
for (i in unique(ref$region)) {
ref[ref$region == i,]$fCO2 <- c(285.78, cumsum(ref[ref$region == i, ]$oi[-1]) + 285.78)
}
#CO2aq
ref$CO2aq <- exp(ref$lnK0)*ref$fCO2
# Laws expression for a given year
ref$laws.current <- with(ref, laws.fun(e1=1, e2=26.5, e.1=1, p=1.5e-5, laws.CO2=CO2aq, sst=sst, r=r, b=0.2))
data$order <- seq(1,length(data$id),1)
data <- merge(data, ref, c("region", "year"), sort = F)
dat1850 <- ref[ref$year==1850,]
dat1850$laws1850 <- with(dat1850, laws.fun(e1=1, e2=26.5, e.1=1, p=1.5e-5, laws.CO2=CO2aq, sst=sst, r=r, b=0.2))
data <- merge(data, dat1850[,c("region", "laws1850")], "region", sort = F)
dat.correct.to <- ref[ref$year==correct.to,]
dat.correct.to$laws.correct.to <- with(dat.correct.to, laws.fun(e1=1, e2=26.5, e.1=1, p=1.5e-5, laws.CO2=CO2aq, sst=sst, r=r, b=0.2))
data <- merge(data, dat.correct.to[,c("region", "laws.correct.to")], "region", sort = F)
data$Laws.cor <- with(data, (laws.current - laws1850) - (laws.correct.to - laws1850))
SuessR.out <- data[,c("id","year", "d13c", "Laws.cor")]
names(SuessR.out)[3] <- "d13c.uncor"
SuessR.out$Suess.cor <- with(data, up.con*exp((year-1850)*0.027)
- up.con*exp((correct.to-1850)*0.027))
SuessR.out$net.cor <- SuessR.out$Suess.cor + SuessR.out$Laws.cor
SuessR.out$d13c.cor <- data$d13c + SuessR.out$net.cor
SuessR.out <- SuessR.out[order(data$order),]
print(cbind(SuessR.out[,1:3], round(SuessR.out[,4:7],3)))
}
# Function for calculating regional uptake constants
#' @encoding UTF-8
#' @title Calculate a regional carbon dioxide uptake constant from empirical DIC data
#' @description The reg.uptake() function calculates the regional uptake constant required to modify the global Suess effect
#' curve to be region-specific.
#' @param year1 The year in which samples were collected for the first set of DIC \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C \cr observations.
#' @param year2 The year in which samples were collected for the second set of DIC \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C \cr observations.
#' @param d13c.change The observed change in DIC \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C between the two years.
#' @details This function calculates the regional uptake constant used to modify the global Suess effect curve to make it
#' specific to a region. This method requires \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C values for DIC from the area of interest from two different time points,
#' ideally separated by at least a decade. The function uses the magnitude (i.e., absolute value) of the observed change in
#' DIC \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C ('d13c.change'), the year of the first observation ('year1'), and the year of the second observation ('year2')
#' to calculate the regional uptake constant. This value can then be supplied as part of the 'custom.region.data' argument
#' (filling a column titled 'up.con') for the SuessR.custom() function.
#' @return Returns a numerical value representing the regional uptake constant. This value can be supplied as part of the
#' 'custom.region.data' (filling a column titled 'up.con') argument for the SuessR.custom() function.
#' @examples
#' year1 <- 1970
#' year2 <- 1980
#' d13c.change <- 0.3
#'
#' reg.uptake(year1 = year1, year2 = year2, d13c.change = d13c.change)
#'
#' @export
reg.uptake <- function(year1, year2, d13c.change) {
reg.up.const <- round(d13c.change/(exp((year2-1850)*0.027) - exp((year1 - 1850)*0.027)),3)
print(reg.up.const)
}
ctclark/SuessR documentation built on March 20, 2021, 5:23 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions reg.uptake SuessR.custom SuessR laws.fun
Documented in reg.uptake SuessR SuessR.custom
# SuessR Package
globalVariables("SuessR.reference.data")
# Function to calculate the Laws correction. Used within the SuessR() and SuessR.custom() functions.
laws.fun <- function(e1, e2, e.1, laws.CO2, sst, r, b, p) {
e2 + e1 - e.1 - (1/(1+laws.CO2*((((p)*(86400)*4*pi*(r*10^-6)^2))/((0.89*1.47^((sst + 273.15 - 303.15)/10))*((0.216*((4/3)*pi*r^3)^0.939)*10^-12)*(1+b)))))*((e2 - e.1)/(b+1))
}
# Function to calculate and apply Suess and Laws corrections.
#' @encoding UTF-8
#' @title Calculate Suess and Laws corrections for stable carbon isotope data from a built-in region
#' @description Generates region-specific Suess, Laws, and net (Suess + Laws) corrections for stable carbon isotope data
#' \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data input by the user. The net correction is then used to calculate the corrected \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data, which are supplied
#' in the output. This function is specifically for data from regions currently built into SuessR
#' ("Bering", "Aleutians","Gulf of Alaska", and "Subpolar North Atlantic", as of Version 0.1.3).
#' \if{html}{\figure{Built-in-regions.png}{options: width=600 alt="Built-in regions"}}
#' \if{latex}{\figure{Built-in-regions.png}{options: width=6in}}
#' @param data A matrix or data frame including columns containing sample ID ('id'), year of sample collection ('year'),
#' region of sample collection ('region'), and uncorrected \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C
#' data ('d13c').
#' @param correct.to The year to which the \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data will be corrected. Defaults to 1850 (Industrial Revolution).
#' @details The SuessR() function uses the year and region of sample collection to calculate Suess, Laws, and net (Suess + Laws)
#' corrections for \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data from marine organisms. The Suess Correction represents the change in \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C values of dissolved
#' inorganic carbon (DIC) in the surface ocean, and is calculated using an exponential function, calibrated to the global decline
#' DIC \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C values. Surface waters in different regions exhibit varying rates of CO\ifelse{html}{\out{<sub>2</sub>}}{\eqn{_2}} uptake from the atmosphere as a result of
#' differing water mass properties and residence time at the surface. Thus, a region-specific modifier has been
#' calculated from empirical observations of changes in DIC \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C values through time and is applied to each region. The Laws
#' correction accounts for changes in stable carbon isotope fractionation during CO\ifelse{html}{\out{<sub>2</sub>}}{\eqn{_2}} uptake and photosynthesis by phytoplankton,
#' and is impacted by changes in aqueous CO\ifelse{html}{\out{<sub>2</sub>}}{\eqn{_2}} concentrations, temperature, and salinity, as well as community growth rates,
#' average cell diameter, average organic carbon content of phytoplankton cells, permeability of the phytoplankton plasmalemma to CO\ifelse{html}{\out{<sub>2</sub>}}{\eqn{_2}},
#' and the ratio of net diffusional loss of CO\ifelse{html}{\out{<sub>2</sub>}}{\eqn{_2}} to carbon fixation. Historic observations and reconstructions of atmospheric CO\ifelse{html}{\out{<sub>2</sub>}}{\eqn{_2}} concentrations,
#' sea surface temperature, and sea surface salinity are used to make these calculations. See references for more details.
#'
#' This function allows users to correct their data to a specific year, using the 'correct.to' argument. This argument defaults
#' to AD1850, which represents onset of the Industrial Revolution and the initiation of the large-scale combustion of fossil
#' fuels that created the Suess Effect. Correcting data to 1850 allows them to be compared to data from any samples collected during
#' or before 1850 (including archaeological samples), as well as to any other samples Suess corrected to the year 1850. Users
#' examining contemporary datasets may wish to correct their data to another year. For example, a user examining a dataset
#' of \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data spanning the years 1970-2010 might choose to correct all their data back to the year 1970, to correct the
#' older samples forward to 2010, or to correct all the samples to 1990, the middle of the time series. In any of these scenarios,
#' corrected data would be comparable to one another. When choosing a value for 'correct.to' consideration should be given to
#' the comparability of the user's data to existing data sets, clarity of presentation of results (i.e., clear statements of the
#' year to which the data were corrected), and reproducibility of results (i.e., presentation of both uncorrected and corrected
#' data so future users can repeat analyses or correct the raw data to another year).
#' @return The output of this function is a data frame that includes the sample ID ('id'), year ('year'),
#' uncorrected \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C ('d13c.uncor'), Laws correction ('Laws.cor'), Suess Correction ('Suess.cor'), net correction
#' ('net.cor'), and corrected \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C ('d13c.cor') for each sample. The corrected \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data is equal to the uncorrected \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data
#' plus the net correction. The units for all values are the standard 'per mil' used for \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data.
#' @examples
#' example.data <- data.frame(id = c("Sample 1", "Sample 2", "Sample 3", "Sample 4",
#' "Sample 5", "Sample 6", "Sample 7", "Sample 8"),
#' year = c(2017, 2017, 2017, 2017, 1977, 1977, 1977, 1977),
#' d13c = c(-12, -12, -12, -12, -12, -12, -12, -12),
#' region = c("Bering Sea", "Aleutian Islands",
#' "Gulf of Alaska", "Subpolar North Atlantic",
#' "Bering Sea", "Aleutian Islands",
#' "Gulf of Alaska", "Subpolar North Atlantic"))
#'
#' SuessR(data = example.data)
#' @references Clark, C.T., M.R. Cape, M.D. Shapley, F.J. Mueter, B.P. Finney, and N. Misarti. (In Prep) SuessR: Regional Suess
#' and Laws corrections for \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C from marine organisms.
#'
#' Clark, C.T., L. Horstmann, A. de Vernal, A.M. Jensen, and N. Misarti. (2019) Pacific walrus diet across 4000 years of
#' changing sea ice conditions. \emph{Quaternary Research}, 1-17.
#'
#' Misarti, N., B. Finney, H. Maschner, and M.J. Wooller. (2009) Changes in northeast Pacific marine ecosystems over the last 4500
#' years: evidence from stable isotope analysis of bone collagen from archaeological middens. \emph{The Holocene}, 19:8. 1139-1151.
#'
#' @importFrom stats na.omit
#' @export
SuessR <- function(data, correct.to = 1850) {
# Make sure the data is a dataframe
if(is.matrix(data)) data <- as.data.frame(data)
# Make input data column names are correct
if(any(!(colnames(data) %in% c("id", "year", "region", "d13c")))) stop(paste("Unrecognized input data:", unique(colnames(data)[!(colnames(data) %in% c("id", "year", "region", "d13c"))]), " "))
# Screen input data to check for (and remove) missing values, with warning
if(anyNA(data)) warning("Rows with missing values were deleted")
data <- na.omit(data)
# Screen input data to make sure all years provided are integers between 1850 and the current year -1
current.year <- as.numeric(substr(date(), 21,24))-1
if(any(!(data$year %in% 1850:current.year))) stop(paste("Year must be an integer between 1850 and ", current.year))
# Screen input data to make sure the 'correct.to' year is an integer between 1850 and the current year -1
current.year <- as.numeric(substr(date(), 21,24))-1
if(any(!(correct.to %in% 1850:current.year))) stop(paste("correct.to must be an integer between 1850 and ", current.year))
# Screen input data to check for non-negative d13c values
if(any(data$d13c >= 0)) warning(paste("Some d13c value(s) are non-negative"))
# Screen input data to check for (and remove) missing values, with warning
if(any(!(data$region %in% SuessR.reference.data$region))) stop(paste("Unrecognized region:", unique(data$region[!(data$region %in% SuessR.reference.data$region)]), " "))
# Create output data frame
SuessR.out <- data.frame(id = data$id, year = data$year)
ref <- SuessR.reference.data
#ln(K0)
ref$lnK0 <- with(ref,
-58.0931 + 90.5069 * (100/(sst + 273.15))
+ (22.294 * log((sst + 273.15)/100))
+ (sss * (0.027766 + -0.025888 * ((sst+273.15)/100)
+ 0.005058 * ((sst + 273.15)/100)^2)))
#Ocean Increase
ref$oi <- NA
for(i in unique(ref$region)) {
sub.i <- ref[ref$region == i,]
n <- nrow(sub.i)
ref[ref$region == i, ]$oi <- c(NA, with(sub.i, (CO2atm[-1] - CO2atm[-n]) * sub.i$Cp[1]))
}
#~f(CO2)ocean
ref$fCO2 <- NA
for (i in unique(ref$region)) {
ref[ref$region == i,]$fCO2 <- c(285.78, cumsum(ref[ref$region == i, ]$oi[-1]) + 285.78)
}
#CO2aq
ref$CO2aq <- exp(ref$lnK0)*ref$fCO2
# Laws expression for a given year
ref$laws.current <- with(ref, laws.fun(e1=1, e2=26.5, e.1=1, p=1.5e-5, laws.CO2=CO2aq, sst=sst, r=r, b=0.2))
data$order <- seq(1,length(data$id),1)
data <- merge(data, ref, c("region", "year"), sort = F)
dat1850 <- ref[ref$year==1850,]
dat1850$laws1850 <- with(dat1850, laws.fun(e1=1, e2=26.5, e.1=1, p=1.5e-5, laws.CO2=CO2aq, sst=sst, r=r, b=0.2))
data <- merge(data, dat1850[,c("region", "laws1850")], "region", sort = F)
dat.correct.to <- ref[ref$year==correct.to,]
dat.correct.to$laws.correct.to <- with(dat.correct.to, laws.fun(e1=1, e2=26.5, e.1=1, p=1.5e-5, laws.CO2=CO2aq, sst=sst, r=r, b=0.2))
data <- merge(data, dat.correct.to[,c("region", "laws.correct.to")], "region", sort = F)
data$Laws.cor <- with(data, (laws.current - laws1850) - (laws.correct.to - laws1850))
SuessR.out <- data[,c("id","year", "d13c", "Laws.cor")]
names(SuessR.out)[3] <- "d13c.uncor"
SuessR.out$Suess.cor <- with(data, up.con*exp((year-1850)*0.027)
- up.con*exp((correct.to-1850)*0.027))
SuessR.out$net.cor <- SuessR.out$Suess.cor + SuessR.out$Laws.cor
SuessR.out$d13c.cor <- data$d13c + SuessR.out$net.cor
SuessR.out <- SuessR.out[order(data$order),]
print(cbind(SuessR.out[,1:3], round(SuessR.out[,4:7],3)))
}
# A function that allows users to add custom regions for the SuessR calculations
#' @encoding UTF-8
#' @title Calculate Suess and Laws corrections for stable carbon isotope data from a custom region
#' @description The SuessR.custom() function generates region-specific Suess, Laws, and net (Suess + Laws) corrections for
#' \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data input by the user. The net correction is then used to calculate the corrected \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data,
#' which are supplied in the output. This function is specifically for data from regions not currently built into SuessR
#' (i.e., different from "Bering", "Aleutians", "Gulf of Alaska", and "Subpolar North Atlantic", as of Version 0.1.3). Or
#' for scenarios in which users want to change the default parameters used by SuessR for the built-in regions (e.g., averge phytoplankton cell radius).
#' @param data A dataframe including sample ID, year of sample collection, uncorrected \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C
#' data, and region. This function is specifically for data from regions not
#' currently built into SuessR. Columns must be named 'id', 'year', 'd13c', and 'region'.
#' @param custom.region.data A data frame containing environmental data for the custom region from which the
#' samples originated. Must contain columns titled 'year', 'region', 'r', 'sst', 'sss', 'CO2atm', 'up.con', and 'Cp'.
#' See details for information on how to supply these parameters appropriately.
#' @param correct.to The year to which the \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data will be corrected. Defaults to 1850, pre-Suess effect.
#' @details The SuessR.custom() allows users to calculate and apply Suess and Laws corrections to \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data from marine
#' organisms collected in a region not currently built into the SuessR package. In Version 0.1.3, the
#' built-in regions are the Bering Sea ('Bering'), the Aleutian archipelago ('Aleutians'), and the Gulf of Alaska ('Gulf of
#' Alaska'). Because the Suess and Laws corrections require region-specific environmental data from 1850-present, users must
#' supply these data using the 'custom.region.data' argument. Once these data have been supplied, this function calculates
#' the corrections exactly like the SuessR() function. See the built-in 'SuessR.reference.data' object for an example
#' template for the data, as well as the references for detailed information on how to compile and supply the appropriate
#' data to the SuessR.custom() function. After compiling these data, please consider sharing them with the package authors
#' (ctclark 'at' alaska.edu) to be included as built-in regions in future version of this package.
#'
#' As with the SuessR() function, SuessR.custom() allows users to correct their data to a specific year, using the 'correct.to'
#' argument. This argument defaults to AD1850, which represents onset of the Industrial Revolution and the initiation of the
#' large-scale combusion of fossil fuels that created the Suess Effect. Correcting data to 1850 allows them to be compared to
#' data from any samples collected during or before 1850 (including archaeological samples), as well as to any other samples
#' Suess corrected to the year 1850. Users examining contemporary datasets may wish to correct their data to another time year.
#' For example, a user examining a dataset of \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data spanning the years 1970-2010 might choose to correct all their data back
#' to the year 1970, to correct the older samples forward to 2010, or to correct all the samples to 1990, the middle of the time
#' series. In any of these scenarios, corrected data would be comparable to one another. When choosing a value for 'correct.to'
#' consideration should be given to the comparability of the user's data to existing data sets, clarity of presentation of
#' results (i.e., clear statements of the year to which the data were corrected), and reproducibility of results (i.e.,
#' presentation of both uncorrected and corrected data so future users can repeat analyses or correct the raw data to another year).
#' @return The output of this function is a data frame that includes the sample ID ('id'), year ('year'),
#' uncorrected \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C ('d13c.uncor'), Laws correction ('Laws.cor'), Suess Correction ('Suess.cor'), net correction
#' ('net.cor'), and corrected \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C ('d13c.cor') for each sample. The corrected \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data is equal to the uncorrected \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data
#' plus the net correction. The units for all values are the standard 'per mil' used for \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C data.
#' @examples
#' example.region.data <- data.frame(year = seq(from = 1850, to = 2020, by = 1),
#' region = rep("Example Region", 171),
#' r = rep(5, 171),
#' sst = seq(5.9, 6.6, 171),
#' sss = seq(32.3, 32.7, 171),
#' CO2atm = SuessR.reference.data$CO2atm,
#' up.con = rep(0.014, 171),
#' Cp = rep(0.41, 171))
#'
#' example.custom.data <- data.frame(id = c("Sample 1", "Sample 2",
#' "Sample 3", "Sample 4"),
#' year = c(1850, 1900, 1950, 2000),
#' d13c = c(-12, -12, -12, -12),
#' region = rep("Example Region", 4))
#'
#' SuessR.custom(data = example.custom.data,
#' custom.region.data = example.region.data,
#' correct.to = 1850)
#'
#' @importFrom stats na.omit
#' @export
SuessR.custom <- function(data, custom.region.data, correct.to = 1850) {
# Make sure the data is a dataframe
if(is.matrix(data)) data <- as.data.frame(data)
# Screen input data to check for (and remove) missing values, with warning
if(anyNA(data)) warning("Rows with missing values were deleted")
data <- na.omit(data)
# Make input data column names are correct
if(any(!(colnames(data) %in% c("id", "year", "region", "d13c")))) stop(paste("Unrecognized input data:", unique(colnames(data)[!(colnames(data) %in% c("id", "year", "region", "d13c"))]), " "))
# Screen input data to make sure all years provided are integers between 1850 and the current year -1
current.year <- as.numeric(substr(date(), 21,24))-1
if(any(!(data$year %in% 1850:current.year))) stop(paste("Year must be an integer between 1850 and ", current.year))
# Screen input data to make sure the 'correct.to' year is an integer between 1850 and the current year -1
current.year <- as.numeric(substr(date(), 21,24))-1
if(any(!(correct.to %in% 1850:current.year))) stop(paste("correct.to must be an integer between 1850 and ", current.year))
# Screen input data to check for non-negative d13c values
if(any(data$d13c >= 0)) warning(paste("Some d13c value(s) are non-negative"))
# Create output data frame
SuessR.out <- data.frame(id = data$id, year = data$year)
ref <- rbind(SuessR.reference.data, custom.region.data)
# Screen input data to check for (and remove) missing values, with warning
if(any(!(data$region %in% ref$region))) stop(paste("Unrecognized region: "), data$region[!(data$region %in% ref$region)])
#ln(K0)
ref$lnK0 <- with(ref,
-58.0931 + 90.5069 * (100/(sst + 273.15))
+ (22.294 * log((sst + 273.15)/100))
+ (sss * (0.027766 + -0.025888 * ((sst+273.15)/100)
+ 0.005058 * ((sst + 273.15)/100)^2)))
#Ocean Increase
ref$oi <- NA
for(i in unique(ref$region)) {
sub.i <- ref[ref$region == i,]
n <- nrow(sub.i)
ref[ref$region == i, ]$oi <- c(NA, with(sub.i, (CO2atm[-1] - CO2atm[-n]) * sub.i$Cp[1]))
}
#~f(CO2)ocean
ref$fCO2 <- NA
for (i in unique(ref$region)) {
ref[ref$region == i,]$fCO2 <- c(285.78, cumsum(ref[ref$region == i, ]$oi[-1]) + 285.78)
}
#CO2aq
ref$CO2aq <- exp(ref$lnK0)*ref$fCO2
# Laws expression for a given year
ref$laws.current <- with(ref, laws.fun(e1=1, e2=26.5, e.1=1, p=1.5e-5, laws.CO2=CO2aq, sst=sst, r=r, b=0.2))
data$order <- seq(1,length(data$id),1)
data <- merge(data, ref, c("region", "year"), sort = F)
dat1850 <- ref[ref$year==1850,]
dat1850$laws1850 <- with(dat1850, laws.fun(e1=1, e2=26.5, e.1=1, p=1.5e-5, laws.CO2=CO2aq, sst=sst, r=r, b=0.2))
data <- merge(data, dat1850[,c("region", "laws1850")], "region", sort = F)
dat.correct.to <- ref[ref$year==correct.to,]
dat.correct.to$laws.correct.to <- with(dat.correct.to, laws.fun(e1=1, e2=26.5, e.1=1, p=1.5e-5, laws.CO2=CO2aq, sst=sst, r=r, b=0.2))
data <- merge(data, dat.correct.to[,c("region", "laws.correct.to")], "region", sort = F)
data$Laws.cor <- with(data, (laws.current - laws1850) - (laws.correct.to - laws1850))
SuessR.out <- data[,c("id","year", "d13c", "Laws.cor")]
names(SuessR.out)[3] <- "d13c.uncor"
SuessR.out$Suess.cor <- with(data, up.con*exp((year-1850)*0.027)
- up.con*exp((correct.to-1850)*0.027))
SuessR.out$net.cor <- SuessR.out$Suess.cor + SuessR.out$Laws.cor
SuessR.out$d13c.cor <- data$d13c + SuessR.out$net.cor
SuessR.out <- SuessR.out[order(data$order),]
print(cbind(SuessR.out[,1:3], round(SuessR.out[,4:7],3)))
}
# Function for calculating regional uptake constants
#' @encoding UTF-8
#' @title Calculate a regional carbon dioxide uptake constant from empirical DIC data
#' @description The reg.uptake() function calculates the regional uptake constant required to modify the global Suess effect
#' curve to be region-specific.
#' @param year1 The year in which samples were collected for the first set of DIC \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C \cr observations.
#' @param year2 The year in which samples were collected for the second set of DIC \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C \cr observations.
#' @param d13c.change The observed change in DIC \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C between the two years.
#' @details This function calculates the regional uptake constant used to modify the global Suess effect curve to make it
#' specific to a region. This method requires \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C values for DIC from the area of interest from two different time points,
#' ideally separated by at least a decade. The function uses the magnitude (i.e., absolute value) of the observed change in
#' DIC \ifelse{html}{\out{δ<sup>13</sup>}}{\eqn{{\delta}^{13}}}C ('d13c.change'), the year of the first observation ('year1'), and the year of the second observation ('year2')
#' to calculate the regional uptake constant. This value can then be supplied as part of the 'custom.region.data' argument
#' (filling a column titled 'up.con') for the SuessR.custom() function.
#' @return Returns a numerical value representing the regional uptake constant. This value can be supplied as part of the
#' 'custom.region.data' (filling a column titled 'up.con') argument for the SuessR.custom() function.
#' @examples
#' year1 <- 1970
#' year2 <- 1980
#' d13c.change <- 0.3
#'
#' reg.uptake(year1 = year1, year2 = year2, d13c.change = d13c.change)
#'
#' @export
reg.uptake <- function(year1, year2, d13c.change) {
reg.up.const <- round(d13c.change/(exp((year2-1850)*0.027) - exp((year1 - 1850)*0.027)),3)
print(reg.up.const)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.