SuessR: Regional corrections for the effects of anthropogenic CO2 on δ13C data from marine organisms

Anthropogenic CO2 emissions associated with fossil fuel combustion have caused declines in baseline oceanic δ13C values. This phenomenon, called the Suess effect, can confound comparisons of marine δ13C data from different years. The Suess effect can be corrected for mathematically; however, a variety of disparate techniques are currently used, often resulting in corrections that differ substantially. SuessR is a free, user‐friendly tool that allows researchers to calculate and apply regional Suess corrections to δ13C data from marine systems using a unified approach. SuessR updates existing methods for calculating region‐specific Suess corrections for samples collected from 1850 to 2020. It also estimates changes in phytoplankton 13C fractionation associated with increasing water temperature and aqueous CO2 concentrations, referred to here as the Laws effect. SuessR version 0.1.3 contains four built‐in regions, including three in the subpolar North Pacific (Bering Sea, Aleutian Islands and Gulf of Alaska) and one North Atlantic region (Subpolar North Atlantic). Users can also supply environmental data for regions not currently built into SuessR to generate their own custom corrections. In 2020, net corrections (Suess + Laws corrections) were as follows—Bering Sea: 1.29‰; Aleutian Islands: 1.30‰, Gulf of Alaska: 1.30‰; and Subpolar North Atlantic: 1.31‰ (compared to a global atmospheric CO2 change of ~2.43‰ across the same period). For samples collected in 2020, the net correction exceeds instrumental error (±0.2‰) when making comparisons across only eight years (i.e. 2013–2020). The magnitude of the Suess effect calculated by SuessR aligns with published estimates, whereas the Laws effect is smaller than previously calculated, resulting from updated estimates of average community cell sizes, growth rates and permeability of phytoplankton plasmalemmas (the plasma membrane which bounds the cell) to CO2. The increasing magnitude of the Suess effect means this phenomenon is no longer only of concern to historical ecologists, but now affects contemporary ecological studies using δ13C data. This highlights the importance of a unified approach for generating Suess corrections. The SuessR package provides a customizable tool that is simple to use and will improve the interpretability and comparability of future stable isotopic studies of marine ecology.


| INTRODUC TI ON
Anthropogenic emissions of carbon dioxide, primarily through combustion of fossil fuels, have altered the baseline δ 13 C values of the planet's atmosphere and oceans since the Industrial Revolution (Keeling, 1979). This δ 13 C decline is commonly referred to as the Suess effect, and has been increasing in magnitude exponentially, in tandem with the exponential increase in 13 C-depleted CO 2 emissions (Bacastow et al., 1996). The shifting δ 13 C baseline is incorporated into biological systems by way of photosynthetic organisms, which fix ambient CO 2 . These effects propagate throughout food webs, as the stable carbon isotopes of primary producers are incorporated into the tissues of consumers. For ecologists seeking to examine changes in δ 13 C values over time to understand animal migratory movements, food web structure and trophic interactions, as well as animal and plant physiology, the Suess effect has the potential to confound analyses (Misarti et al., 2009).
In addition to the shifts in δ 13 C that characterize the Suess effect, increasing atmospheric CO 2 concentrations affect plant physiology in ways that may impact stable carbon isotope fractionation during photosynthesis (Keeling et al., 2017). In marine and aquatic systems, 13 C fractionation by phytoplankton changes in relation to water temperature and aqueous CO 2 concentrations (Laws et al., 2002). Thus, warming of the global oceans observed in recent decades, as well as substantial increases in atmospheric and aqueous CO 2 concentrations during this same period, has also altered baseline δ 13 C values in marine ecosystems. The impacts of water temperature (via community composition and growth rates) and aqueous CO 2 concentrations on stable carbon isotope fractionation by phytoplankton were described mathematically by Laws et al. (2002), and are hereafter referred to as the Laws effect.
Efforts have been made to quantify and mathematically correct for the Suess and Laws effects, primarily by historical ecologists, palaeoecologists and archaeologists comparing δ 13 C values across broad timeframes. These researchers typically correct historic and modern values back to preindustrial times (i.e. the year 1850, the onset of the Industrial era and the beginning of exponential increases in anthropogenic CO 2 emissions into the atmosphere and oceans; Ruddiman, 2013). However, as the burning of fossil fuels has continued to produce exponentially increasing releases of anthropogenic CO 2 (Bacastow et al., 1996), the year-to-year changes in the magnitude of the Suess effect (and to a lesser extent, the Laws effect) have increased substantially. As of 2020, the summed magnitude of the Suess and Laws corrections for δ 13 C values from the Bering Sea, calculated using the methods described in this paper, exceeds the typically reported value for instrumental precision for isotope ratio mass spectrometers (±0.2‰) when comparing samples across just 8 years (i.e. comparing samples from 2020 to samples from 2013; Figure 1). Thus, the Suess and Laws effects introduce detectable error in comparisons of modern stable carbon isotope data collected over less than a decade. This suggests that corrections for the Suess and Laws effects are no longer only the concern of historical ecologists, palaeoecologists and archaeologists, and that these corrections should be incorporated into contemporary ecological studies comparing isotope values over spans of 8 years or more. Although calculation of mathematical Suess corrections has become relatively commonplace, the methods used vary widely and can produce substantially different results ( Figure 1), highlighting the need for a unified approach to calculate Suess and Laws corrections for δ 13 C data.
The SueSSR package is a tool for calculating and applying mathematical corrections for the Suess and Laws effects. These corrections are region-specific, accounting for the complex spatial variability in seawater circulation, surface residence time, water temperature and biological production that impact CO 2 uptake by the oceans (Eide et al., 2017), as compared with the relatively homogeneous and well-mixed global atmosphere. The purpose of this R package is to provide a simple and effective tool that allows ecologists to apply these corrections to δ 13 C data from marine systems using a single, unified methodology. Additionally, users will be able to supply their own region-specific parameters to the Suess and Laws corrections, and are encouraged to submit these data for incorporation into the rates and permeability of phytoplankton plasmalemmas (the plasma membrane which bounds the cell) to CO 2 .
4. The increasing magnitude of the Suess effect means this phenomenon is no longer only of concern to historical ecologists, but now affects contemporary ecological studies using δ 13 C data. This highlights the importance of a unified approach for generating Suess corrections. The SueSSR package provides a customizable tool that is simple to use and will improve the interpretability and comparability of future stable isotopic studies of marine ecology.

K E Y W O R D S
carbon dioxide, carbon isotopes, fractionation, phytoplankton, Suess correction,

Suess effect
SueSSR package, thereby increasing the number of built-in regions available in future versions.

| The Suess correction
Marine ecologists have used a variety of approaches to mathematically correct for the Suess effect. The magnitude of these corrections can vary widely, and a single approach often gains popularity within an individual sub-field (e.g. seabird research), making it more difficult to compare the results of studies performed on different taxa. The most basic Suess correction involves adding a fixed value (e.g. 1‰) to the δ 13 C values of modern samples to make them comparable to archaeological or historic specimens that lived prior to the Industrial Revolution (e.g. Conrad et al., 2018;Elliott Smith et al., 2020;Halffman et al., 2015;Newsome et al., 2004;Zangrando et al., 2016). This method of correction is simple and straightforward, but requires that modern samples represent a 'snapshot', rather than a time series, and that they be collected within a few years of one another. Furthermore, it requires a good estimate of the magnitude of the Suess effect in the year(s) of sample collection. The use of a 1‰ correction by papers published as early as 2004 and as late as 2016 indicates that there is room for improvement in these estimates. To allow for comparison of samples collected at different times during the Industrial era, researchers have used linear estimates of the rate of change of δ 13 C values associated with the Suess effect, typically based on δ 13 C values of dissolved inorganic carbon (DIC) in a particular region or ocean basin (e.g. -0.018‰ per year across the entire North Atlantic basin; Quay et al., 2007), to calculate annual corrections (e.g. Espinasse et al., 2019;Matthews & Ferguson, 2014;Ramos et al., 2020;Rossman et al., 2013;Soto et al., 2018). An improved version of this approach incorporates an increase in the rate of δ 13 C decline after a specific year (e.g. 1950), to better approximate the exponential shape of the Suess effect curve (e.g. Alter et al., 2012;Bond & Lavers, 2014;Farmer & Leonard, 2011;Sun et al., 2019;Vokhshoori et al., 2019). Both of these methods allow for region-specific Suess corrections; however, there are also added complications. Researchers will typically correct all samples back to the first year of their time series (e.g. 1930) or forward to the last year (e.g. 2010), thereby generating the Suess corrected δ 13 C values that are comparable to one another, but not to pre-Industrial samples or to results produced from other studies that are corrected to a different year. When used to correct samples back to preindustrial δ 13 C levels (i.e. back to ~1850), these approaches tend to substantially overestimate the magnitude of the Suess effect (Figure 1), though they may produce more realistic corrections when applied to a restricted period of time.
Other Suess corrections have been generated using nonlinear models based on measured changes in environmental δ 13 C values.
Researchers working in lakes have used higher-order polynomial models based on atmospheric trends to estimate the magnitude of the atmospheric Suess effect, assuming these systems are close to equilibrium with the atmosphere (e.g. Verburg, 2007). Some studies have applied these atmospheric corrections to marine systems by applying a lag to account for the incorporation of atmospheric CO 2 into surface waters (e.g. Vales et al., 2020), which is likely a better approximation of the Suess effect than provided by linear models; however, it is unclear whether the 10-year lag appropriately captures the disequilibrium between the atmospheric and oceanic Suess effect (Figure 1). Hilton et al. (2006) used values from the published literature to develop a more complex Suess correction, which accounts for the observed exponential decline in δ 13 C values of dissolved inorganic carbon in the global oceans associated with the absorption of anthropogenically produced CO 2 . This correction also included a constant accounting for regional differences in CO 2 uptake by the ocean. Misarti et al. (2009) subsequently updated and adapted this correction for the Gulf of Alaska, and corrected an error in the original equation published by Hilton et al. (2006), where an addition symbol had been substituted for a multiplication symbol. The resulting equation (Misarti et al., 2009) is: In this equation, a is a constant reflecting the maximum observed rate of δ 13 C decline in surface waters DIC in a specific region, b is the year of sample collection minus 1850 (the onset of the Industrial Revolution) and 0.027 is the parameter value obtained by Hilton et al. (2006) after fitting an exponential curve to the global ocean δ 13 C data from 1945 to 1997 published by Gruber et al. (1999).
The resulting Suess effect correction (a negative value) may then be subtracted from the uncorrected δ 13 C data to obtain the Suesscorrected δ 13 C value. This equation has since been used by numerous studies to calculate regional Suess corrections (e.g. Clark et al., 2019;Guiry et al., 2020;Harris et al., 2020;Kochi et al., 2018).

| The Laws correction
The second correction calculated and applied by SuessR was developed by Laws et al. (2002), and takes into account shifts in stable carbon isotope fractionation by phytoplankton associated with changes in water temperature, aqueous CO 2 concentrations and phytoplankton physiology. The equation for estimating the difference between aqueous CO 2 δ 13 C values and those of the products of photosynthesis by phytoplankton (i.e. stable carbon isotope fractionation by phytoplankton, ε p ) is: In this equation, ε 1 represents the isotopic fractionation associated with diffusion of DIC into the cell, whereas ε −1 represents the isotopic fractionation associated with diffusion of DIC out of the cell. Both of these terms are assumed to be constant at 1‰, and to cancel one another out (Laws et al., 2002); thus, their appearance at the beginning of the above equation serves primarily as a reminder of this assumption.
The isotopic fractionation associated with carboxylation, represented in this equation by ε 2 , is assumed to be constant at 26.5‰, the median of the values presented by Laws et al. (2002). β is the ratio of net diffusional loss of CO 2 to carbon fixation, and is assumed to be constant. P reflects the permeability of the plasmalemma (the plasma membrane which bounds the cell) to CO 2 (p) multiplied by the surface area of the cell (thus, P = p × cell surface area, with p in m/day), C the organic carbon content of the cell (in gC/cell) and μ is the average growth rate (in day −1 ) of phytoplankton in a given region and year. Both P and C are proportional to cell size, with P scaling to a cell's surface area and C to cell volume. [CO 2 ] aq (in µmol/kg) represents the concentration of aqueous CO 2 in the seawater.
Estimates of [CO 2 ] aq are calculated by multiplying the estimated fugacity of the CO 2 in the ocean (fCO 2 ocean , in µatm) by the solubility of CO 2 in seawater at a given temperature and salinity (K 0 ). Calculating annual estimates of fCO 2 ocean first requires comparing the regional rates of change of fCO 2 atmosphere and fCO 2 ocean to determine the proportion of atmospheric CO 2 concentration increase exhibited by the ocean in a given year. For example, if fCO 2 atmosphere in a particular region increases by 10 µatm and fCO 2 ocean by 6 µatm over the course of a decade, the annual increase in surface seawater fugacity (ΔfCO 2 ocean ) for this region can be estimated by multiplying measured changes in fCO 2 atmosphere ΔfCO 2 atmosphere by 0.6. We refer to this number as the proportional rate of change (C P ). Annual estimates of fCO 2 ocean can thus be calculated using the following equation: This approach requires an assumed starting value for fCO 2 ocean at the beginning of the record (fCO 2 ocean t=0 ). Additionally, because records of atmospheric CO 2 concentrations are presented in parts per million (ppm), this approach assumes that CO 2 behaves as an ideal gas and that fCO 2 atmosphere is directly proportional to atmospheric CO 2 concentrations.
Annual estimates of K 0 can be generated using the following equation from Weiss (1974): In this equation, absolute temperature (T, in degrees Kelvin) and sea surface salinity (SSS), along with constants A 1-3 and B 1-3 , are used to calculate the natural log of K 0 in ln(moles/l·atm). Once these terms have been calculated, they can be used to generate yearly estimates of [CO 2 ] aq using the following equation: To calculate the Laws correction, which is the overall change in ε p by phytoplankton since 1850, the estimated stable carbon isotope fractionation from 1850 is subtracted from that of the year in which a sample was collected: Sea surface temperature (SST) and aqueous CO 2 concentrations have opposite effects on ε p , with increases in SST leading to smaller ε p values, and increases in [CO 2 ] aq resulting in higher values for ε p . Since 1850, increases in aqueous CO 2 concentrations have had more of an effect on ε p than have increasing water temperatures, and ε p has increased across the last ~170 years. Thus, the Laws correction is typically a positive number added to the sample's δ 13 C value to account for this increase in ε p (because increasing 13 C fractionation through time means more negative δ 13 C values in later years) and provide the Laws-corrected data.

| Overview of SueSSR
Version 0.1.3 of SuessR includes four built-in regions for which users may calculate Suess and Laws corrections and apply them to their data. These include three North Pacific regions (Bering Sea, Aleutian Islands, Gulf of Alaska; Figure 2; Table S1) and one North Atlantic region (Subpolar North Atlantic; Figure 2; Table S1). Additionally, the package includes the SuessR.custom() function, which allows users to supply their own regional parameters for the Suess and Laws Corrections. It is our hope that users will submit these data to the authors of the package for inclusion in future versions of SuessR.
The updates to the equation used to generate the Laws correction and the introduction of the function for calculating the regional uptake constant for the Suess correction were primarily aimed at facilitating the process of using the SueSSR package for custom regions.
Users choosing to parameterize the Suess and Laws corrections for custom regions should have a thorough understanding of the physical, chemical and biological systems in the region of interest, and careful consideration is required to avoid erroneous results. Users intending to do so must make a number of important considerations pertaining to their region(s) of interest and the organisms sampled by their study. SuessR is available on CRAN (https://cran.r-proje ct.org/ web/packa ges/Suess R/index.html) and can be installed command 'install.packages('SuessR')' in the R console. A web-based Shiny app that allows users to interact with SuessR using a graphical user interface is also available (https://suessr.shiny apps.io/Suess R/).

| Calculation of the Suess correction
The SueSSR package uses the equation developed by Hilton et al. (2006) and modified by Misarti et al. (2009), but introduces a different method for calculating the regional uptake constant a. To calculate this constant, the observed decline in DIC δ 13 C in a specific region over a given period is typically treated as a linear change and is divided by the number of years over which the change was observed, thus producing an annual rate of δ 13 C change. Taking the linear slope of an exponential curve will, however, produce substantially different results depending on which portion of the curve is used for the calculation. Thus, data from the same region will provide different rates of δ 13 C change if data are examined from 1970 to 1979 and from 2000 to 2009, with the latter producing a larger annual change. To address this problem, the SueSSR package contains a function that derives the regional uptake constant using a different approach. Because the observed change in δ 13 C DIC across a span of time closely approximates the magnitude of the Suess effect across that period, we can solve for the regional uptake constant using the Suess correction equation as follows: In this way, the regional uptake constant for any body of water can be calculated using observed changes in δ 13 C DIC from any time period. This approach assumes that the regional uptake constant for a given body of water has remained the same since 1850. Once this term has been estimated, calculating and applying the Suess correction is straightforward.
The regional uptake constants used for SuessR's built-in regions include an update to those used by Misarti et al. (2009) Table S1 1993 and 2003 using an isopycnal multiple linear regression, which matched well with other estimates from the North Atlantic. This estimate yields a regional uptake constant of -0.013. The identical values calculated for the subpolar North Pacific and subpolar North Atlantic seem reasonable, given that the majority of the variation in the Suess effect in the surface ocean is associated with latitude (Eide et al., 2017;Quay et al., 2003Quay et al., , 2007, although there are some indications that the Suess effect may be slightly more pronounced in the subpolar North Atlantic (Eide et al., 2017).

| Calculation of the Laws correction
Using reconstructions and records of annual atmospheric CO 2 concentrations, sea surface temperature (SST) and sea surface salinity (SSS), yearly estimates of regional aqueous CO 2 concentrations can be generated back to 1850. These data were compiled from various gridded global datasets, using version 0.4.8 of the ReRddapXtRacto package in R (Mendelssohn, 2020). SST data were taken from the  Pacific: 1970Pacific: , 1973Pacific: -1980Pacific: , 1982Pacific: -1983Pacific: , 1985Pacific: -2019Atlantic: 1981Atlantic: -1982Atlantic: , 1989 and 1.00. To be conservative, the calculated value (C p = 1.00) was used for this region, which matches with reports that CO 2 concentrations in North Atlantic surface waters are increasing at about the same rate as the atmosphere (Eide et al., 2017). It is possible that apparent changes in fCO 2 ocean in other regions may also be confounded by changes in water mass properties, and this possibility should be considered when supplying data for custom regions to SuessR. Such changes are typically reported and discussed in the literature, and often manifest as differences between modelled and observed values for rates of δ 13 C DIC change in surface waters.
The SueSSR package uses the equation for calculating p from Laws et al. (2002) as the framework for estimating the Laws correction; however, some components have been updated or amended to make using this equation simpler and to better approximate phytoplankton communities, as opposed to individual species. One of these updated pieces is a different approach for calculating the organic carbon content of the cell, C, than that originally referenced by Laws et al. (2002). This calculation, from Menden-Deuer and Lessard (2000), estimates C based on cell radius, r. The incorporation of this new correction means that only the average cell radius, r (in µ), needs to be supplied to SuessR, simplifying the process for users to generate calculations using data appropriate to their region of interest. To most closely approximate the organic carbon content of the phytoplankton community, SuessR uses the equation This assumes a taxonomically diverse group of protist phytoplankton, but excludes large diatoms with cell volumes >3,000 µm 3 .
Researchers specifically studying diatoms, dinoflagellates or food webs based primarily on these groups might prefer to use the relationship between organic carbon content and cell size specific to these groups, as reported by Menden-Deuer and Lessard (2000); however, SuessR currently does not support the use of these equations, and the Laws correction would need to be calculated manually.
For the majority of phytoplankton communities, however, the equation employed by SuessR should adequately approximate C. Hilton et al. (2006) and Misarti et al. (2009) (Table S2), and the others were held constant at their means. The resulting values were then supplied as parameters to the Laws equation to calculate an estimate of phytoplankton 13 C fractionation (ε p ). This process was repeated 10,000 times for each cell radius (1 µ, 5 µ and 10 µ). The ε p values that resulted from varying each parameter were then plotted (Figure 3). For all parameters, the magnitude of the observed differences increased with increasing cell size. At a radius of 1 µ, none of the parameters produced substantial differences in the overall ε p when varied. At 10 µ, these differences were much larger, and if allowed to vary together, would be expected to compound one another. Finally, the analysis was repeated with all parameters allowed to vary independently.
The results of the sensitivity analysis are important to the interpretation and application of the Laws correction. Varying [CO 2 ] aq and SST resulted in moderate changes to ε p , whereas varying β had little effect on this value (Figure 3). The effects of [CO 2 ] aq and SST were expected, as they were the primary reason for calculating ε p .
Increasing [CO 2 ] aq resulted in larger ε p values, whereas increasing SST resulted in smaller values for ε p . The sensitivity analysis also indicated that ε p is strongly impacted by changes in cell size, which not only changes the calculated value of ε p (larger cells = lower ε p values) but also compounds the effects of varying the other parameters.  Table S2 cell sizes is most apparent in the scenario in which all parameters are allowed to vary, with small cells exhibiting almost no variability in ε p , and large cells varying much more widely. Thus, selecting an average cell radius that best characterizes the phytoplankton communities responsible for the bulk of the primary production in the region of interest is the most important consideration when selecting variables to estimate Laws corrections for a new region. Although published estimates of size distributions within phytoplankton communities are becoming more common in recent years (e.g. Acevedo-Trejos et al., 2014;Bolaños et al., 2020;Laney & Sosik, 2014;Marañón, 2015;Polovina & Woodworth, 2012), there is currently a dearth of empirical data for estimating average cell radius, and more research is needed to better characterize phytoplankton cell size distributions.
The term describing the permeability of the plasmalemma to CO 2 , p, adds additional complications to calculating the Laws correction. The overall carbon fractionation, ε p , is highly sensitive to changes in p ( Figure 3). As p approaches zero, the calculated ε p val-  (Table S3). This value is on the cusp of the sharp declines in ε p visible in Figure 3 that resulted from small decreases in the value of p. The fact that it is slightly larger than the value used by Hilton et al. (2006) and Misarti et al. (2009) is responsible, in part, for the smaller Laws corrections calculated by SuessR, as compared to these two papers.
The value used by SuessR is calculated using data from a variety of phytoplankton groups and is likely a good representation of a community average value for p. That said, the sensitivity of the Laws correction to p, as well as the relative lack of understanding of the natural variability in this term, suggests that it warrants further research. Many of the smaller p values in the published literature (Table S3) would result in unrealistic estimates of ε p , and models suggest that values of p substantially smaller than 1 × 10 -4 m/s are unlikely to be representative of marine phytoplankton species that rely primarily on diffusive uptake of CO 2 (Rau et al., 1996). Future research that generates a better understanding of p will help to constrain this term, and may provide a more realistic estimate of the average p across phytoplankton communities for use in calculating the Laws correction.

| Example application
The  When publishing δ 13 C data corrected by SuessR, it will be important to make clear the year to which the data have been corrected.
SuessR contains an additional function, SuessR.custom(), which allows users to supply their own regional parameters for the Suess and Laws corrections. These supplied parameters are appended to the built-in SuessR.reference.data (Supplementary Dataset SD1) and must be supplied in the same format. To input data for a custom region, the user must first create a data frame with the regional parameters for the Suess and Laws corrections. This data frame can then be called using the 'custom.region.data' argument in the SuessR.custom() function. The name of the custom region should be specified in the 'region' column of the data frame, and this region must be specified when δ 13 C values for this region are supplied for correction. Because the parameters chosen will determine the final values for the Suess and Laws corrections, it is critical that they are carefully selected and vetted before being supplied to SuessR. It is the authors' hope that users will submit the parameters for their custom regions along with the associated references to be incorporated as built-in regions in future versions of SuessR. These custom data can be submitted to the corresponding author of this paper or to the maintainer of the SueSSR package, as laid out in the package documentation.
The reg.uptake() function is used to calculate the regional uptake constant for a region, which can be supplied to the SuessR.custom() function as part of the 'custom.region.data' object. This function requires three arguments: 'year1', 'year2' and 'd13c.change'. The term 'year1' indicates the calendar year in which the first δ 13 C DIC observation was made, whereas 'year2' is the year in which the last δ 13 C DIC measurement was taken. 'd13c.change' represents the magnitude of change in δ 13 C DIC (in ‰) observed between those years. Because of the potential biases and problems associated with snapshot estimates of δ 13 C DIC , many of the more current studies in the published literature use a multiple linear regression approach to estimating the change in δ 13 C DIC through time (Quay et al., 2007). This method more effectively separates natural and anthropogenic δ 13 C DIC changes by accounting for seasonal variability in δ 13 C DIC , and predicts δ 13 C DIC using measured hydrographic data. The regressions typically provide a linear approximation of δ 13 C DIC changes (e.g. −0.018‰ per year from 1990 to 2000). In this example, 'year1' would be 1990, 'year2' would be 2000 and d13c.change would be −0.198‰ (−0.018‰ year −1 × 11 year). The resulting regional uptake constant would be a = −0.015.

| CON CLUS IONS
The SueSSR package provides a widely accessible, customizable and easy to use tool that allows ecologists to apply corrections to account for the influence of anthropogenic CO 2 on marine δ 13 C data. In contrast, the output of the updated Laws equation implemented by SuessR indicates that temperature-and CO 2 -related changes in fractionation by marine phytoplankton are less important than previously thought, although this may apply only to the subpolar regions built into the package, and not at lower latitudes (Young et al., 2013).
The results of some studies indicate that increasing temperature and CO 2 concentrations might cause much larger changes in stable carbon isotope fractionation by marine phytoplankton (ε p ) than predicted by SuessR (Cullen et al., 2001;Young et al., 2013), due primarily to the use of different methods of calculating ε p . Further study is merited to determine which approach better characterizes marine phytoplankton communities, and some updates to SuessR may be required in the future if another method of calculation is found to be more accurate. Improvements and updates to future releases of TA B L E 2 Output from the SuessR() function when supplied the example input data from Table 1. Columns include sample id ('id'), year of sample collection ('year'), uncorrected δ 13 C data ('d13c. uncor', in ‰), the magnitude of the estimated Laws correction ('Laws.cor', in ‰), the magnitude of the estimated Suess correction ('Suess.cor', in ‰), the magnitude of the estimated net correction ('net.cor', in ‰) and the corrected δ 13 C data ('d13c.cor', in ‰) SuessR will rely, in part, on better constraints on model parameters, including r (the average community cell radius) and p (the permeability of the plasmalemma to CO 2 ), which are particularly important to improving the Laws correction. Future versions of the Laws correction might be improved using an iterative process to draw cell sizes from a distribution representing a phytoplankton community, rather than summarizing the community to a single mean value.
Similarly, an update to the long-term changes in δ 13 C DIC modelled by Gruber et al. (1999) will be important to verify that the original Suess correction equation published by Hilton et al. (2006) is still an adequate characterization of the magnitude of the oceanic Suess effect. Furthermore, as some areas of the surface ocean become saturated with CO 2 , the dynamics of the oceanic Suess effect will change in ways that are not currently characterized by the Suess correction. As with any data correction approach, it is important to consider the possibility of bias or error within the correction itself, and presenting both the uncorrected and corrected data is strongly encouraged. This allows researchers and subsequent users of the data to assess the impacts of applying the correction, and to explore alternate hypotheses for observed variability in stable carbon isotope values. Despite these uncertainties and ways in which the Suess and Laws corrections may be refined in the future, the corrections produced by SuessR match with reported estimates of the magnitude of the Suess effect in the surface ocean. This tool has the potential to improve future research by providing a unified approach to calculating and applying Suess corrections for marine δ 13 C data, that is easy to access and free to use. Furthermore, submission of compiled regional environmental data by users will allow for more built-in regions to be added in future versions of SuessR, further increasing its ease of use. Custom regional data can be submitted to the corresponding author or to the SuessR maintainer as designated in the package's documentation.

ACK N OWLED G EM ENTS
This work was funded by the University of Alaska Faculty Initiative Fund, and partial funding was provided to the lead author through NOAA Cooperative Agreements NA15OAR4320063 and NA20OAR4320271. The authors declare no conflicts of interest.

PEER R E V I E W
The peer review history for this article is available at https://publo ns.

DATA AVA I L A B I L I T Y S TAT E M E N T
All data used by SuessR are publicly available in online repositories as specified in the text and references. These datasets include ERSST v5 (SST; https://www.ncdc.noaa.gov/data-acces s/marin eocea n-data/ exten ded-recon struc ted-sea-surfa ce-tempe ratur e-ersst -v5), SODA