Calibration Strategies for Detecting Macroscale Patterns in NEON Atmospheric Carbon Isotope Observations

Carbon fluxes in terrestrial ecosystems and their response to environmental change are a major source of uncertainty in the modern carbon cycle. The National Ecological Observatory Network (NEON) presents the opportunity to merge eddy covariance (EC)‐derived fluxes with CO2 isotope ratio measurements to gain insights into carbon cycle processes. Collected continuously and consistently across >40 sites, NEON EC and isotope data facilitate novel integrative analyses. However, currently provisioned atmospheric isotope data are uncalibrated, greatly limiting ability to perform cross‐site analyses. Here, we present two approaches to calibrating NEON CO2 isotope ratios, along with an R package to calibrate NEON data. We find that calibrating CO2 isotopologues independently yields a lower δ13C bias (<0.05‰) and higher precision (<0.40‰) than directly correcting δ13C with linear regression (bias: <0.11‰, precision: 0.42‰), but with slightly higher error and lower precision in calibrated CO2 mole fraction. The magnitude of the corrections to δ13C and CO2 mole fractions vary substantially by site, underscoring the need for users to apply a consistent calibration framework to data in the NEON archive. Post‐calibration data sets show that site mean annual δ13C correlates negatively with precipitation, temperature, and aridity, but positively with elevation. Forested and agricultural ecosystems exhibit larger gradients in CO2 and δ13C than other sites, particularly during the summer and at night. The overview and analysis tools developed here will facilitate cross‐site analysis using NEON data, provide a model for other continental‐scale observational networks, and enable new advances leveraging the isotope ratios of specific carbon fluxes.


Introduction
The terrestrial carbon cycle and its response to increasing CO 2 concentrations, increased temperatures, and altered water availability is one of the least certain aspects of the climate system (e.g., Bonan et al., 2019;Friedlingstein et al., 2014). Terrestrial carbon cycling is fundamentally linked with the terrestrial energy balance and water cycle, as more than half of the water leaving the land surface during evapotranspiration passes through plants as transpiration (Bowen et al., 2019;Good et al., 2015;Schlesinger & Jasechko, 2014). As a result, several approaches have been developed to improve constraints on terrestrial carbon fluxes. Eddy covariance (EC) has emerged as a popular and powerful technique to investigate fluxes of carbon, water, and energy between the land surface and the atmosphere (e.g., Baldocchi et al., 2001;Loescher et al., 2006), but additional measurements and/or assumptions must be made to estimate constituent fluxes (e.g., Reichstein et al., 2005). For example, partitioning net ecosystem exchange (NEE) from EC measurements into ecosystem gross primary productivity and ecosystem respiration components requires a model (Lasslop et al., 2010;Reichstein et al., 2005). Given current limitations of these models, the resulting component fluxes still exhibit substantial uncertainty (Keenan et al., 2019;Wehr et al., 2016). When paired with EC measurements, complementary measurements of additional tracers that better capture underlying processes driving these component fluxes present a promising approach to reduce these uncertainties.
Stable carbon isotope ratios provide synergistic process-level information on the exchange of carbon between ecosystems and the atmosphere that cannot be gleaned from net carbon fluxes alone. Carbon isotope ratios in atmospheric CO 2 (hereafter δ 13 C-CO 2 ) have been used in a wide array of applications, including: (1) constraining anthropogenic CO 2 emissions, as CO 2 from combusted fossil fuels are comparatively depleted in heavy isotopologues (Pazdur et al., 2007;Sonnerup et al., 1999;Suess, 1955), (2) examining patterns of CO 2 exchange between the land surface, ocean, and the atmosphere, as oceanic fluxes of CO 2 have higher δ 13 C values than terrestrial fluxes (Ciais et al., 1995;Trolier et al., 1996), (3) partitioning NEE into assimilatory and respiratory components (Bowling et al., 2001;Knohl & Buchmann, 2005;Wehr & Saleska, 2015), (4) probing the abundance of C 3 and C 4 photosynthesis and their relative roles in the carbon cycle (Buchmann & Ehleringer, 1998;Ehleringer et al., 1997;Still et al., 2003), (5) examining recycling of respired CO 2 within vegetation canopies (Yakir & Sternberg, 2000), and (6) estimating plant water use efficiency (Baldocchi & Bowling, 2003;Farquhar & Richards, 1984;Medlyn et al., 2017;Seibt et al., 2008). Initial studies of δ 13 C-CO 2 used flask samples (Ciais et al., 1995;Keeling, 1958Keeling, , 1961Trolier et al., 1996), which limited both the spatial and temporal frequency of sampling and our ability to infer process controls and flux magnitudes in terrestrial systems. Advances in laser spectroscopy have allowed for higher frequency measurements and the resolution of δ 13 C-CO 2 atmospheric profiles (Bowling et al., 2005;Griffis et al., 2007;Raczka et al., 2017). These efforts have typically been driven by individual research groups with site-specific protocols and instrumentation; unfortunately, this lack of coordination introduces challenges and uncertainties into cross-site comparison efforts. Furthermore, record lengths vary substantially, and the lack of concurrent multi-year records limits the comparability and information content of existing δ 13 C-CO 2 data sets.
The National Ecological Observatory Network (NEON; Thorpe et al., 2016), supported by the U.S. National Science Foundation, has established a new data stream documenting atmospheric δ 13 C-CO 2 across the United States using standardized instrumentation, data collection and processing protocols. These observations are collocated with EC measurements, as well as periodic measurements of δ 13 C values in other ecosystem pools, such as soil, litter, and plant matter. Measurements of δ 13 C-CO 2 are made at multiple heights on NEON EC towers as part of the EC storage exchange system. At present, NEON calibrates isotopic analyzers once annually, but these observations are otherwise uncalibrated to remove instrumental drift and variation in the base calibration. Yet, three reference standards with known δ 13 C and CO 2 mole fractions are run daily at each analyzer. By making use of these standards, the accuracy of the NEON raw δ 13 C-CO 2 data streams can be corrected to remove this bias associated with drift and tie the absolute values at the network of sites to a common scale. Without this step, robust cross-site data comparison and analyses cannot be guaranteed since the bias may equal or exceed the size of the ecological signal of interest and give rise to spatial difference between sites. In this study, we present calibration strategies and an R package for NEON's δ 13 C-CO 2 measurements and highlight trends apparent across the first 2-3 years of observatory data.
FIORELLA ET AL.

NEON Site Locations and Data Preprocessing
NEON features 20 ecoclimatic domains that span major variations in ecosystem type and climate conditions across the United States. Each domain has a "core" terrestrial site, which is expected to operate for NEON's 30-year mission. In addition, each domain (with the current exceptions of domain 12-Northern Rockies, domain 15-Great Basin, and domain 20-Pacific Tropical) possesses at least one "relocatable" terrestrial site that could potentially move. Both core and relocatable terrestrial sites feature a central instrumented tower, supporting a suite of standard meteorological measurements, an EC system at the tower top, and additional measurements of CO 2 and H 2 O concentrations and isotope ratios at four to eight vertical levels, depending on ecosystem structure. Atmospheric isotope data and EC fluxes are bundled together within NEON data product DP4.00200.001. Additional details on the design, construction, and operation of NEON's EC storage exchange system are provided in the NEON document library .
δ 13 C-CO 2 values are measured at each terrestrial NEON site using a Picarro G2131-i. Ambient air is sampled from each tower intake sequentially, with each height being measured at for 10 min at a time. For each 10 min measurement period, the first minute of each sampling period is excluded to minimize instrumental memory affects associated with switching between intake heights. The Picarro G2131-i analyzer reports δ 13 C values and CO 2 mixing ratios at ∼1 Hz; these measurements are aggregated to 9 min averages during NEON's data processing using the eddy4R software (Metzger et al., 2017;Xu et al., 2019). The return interval for measurements at each height varies by site with tower structure, and ranges between 40 and 80 min under normal operating conditions. As part of the data aggregation from native sensor resolution to 9 min averages, eddy4R calculates quality flags and quality metrics along with data. However, the current iteration of the data flagging framework is overly restrictive for the EC storage exchange system in which δ 13 C-CO 2 measurements are made. These quality flags are currently being reevaluated by NEON, and as a result, NEON quality flags were not utilized in this study. Once per year, the analyzers are shipped from the field to NEON headquarters in Boulder, CO, where the instruments are calibrated using a common laboratory setup. While this approach ensures that measured δ 13 C values throughout the year are generally close to calibrated values, instrumental drift can cause measured values to diverge from the Vienna Pee Dee Belemnite (VPDB) scale throughout the year. Three reference gases with known CO 2 concentrations and δ 13 C values are measured once per day for 10 min each at each NEON site. At present, these measurements are not used to correct measured CO 2 and δ 13 C values for instrumental drift as part of NEON's data processing pipeline. Therefore, the quality and utility of NEON-provisioned δ 13 C-CO 2 measurements could be improved by making use of these daily reference gas measurements.

Carbon Fractionation Background
δ 13 C-CO 2 values provide constraints on carbon cycle processes from ecosystem-to-global spatial scales. At the ecosystem level, variations in δ 13 C-CO 2 primarily reflect the dynamics of net carbon exchange between the ecosystem and the atmosphere via photosynthesis and respiration. These variations in 13 C/ 12 C manifest in part-per-thousand magnitudes, and therefore, isotope ratios are commonly reported in "delta" notation: where R is the molar ratio of 13 CO 2 and 12 CO 2 isotopologues and R VPDB is the heavy-to-light isotope ratio in VPDB (R VPDB = 1.11797 × 10 −2 ) (Griffis et al., 2004;Tans et al., 2017). Carbon isotope ratios for plants are also commonly expressed in terms of an isotopic discrimination for 13 C with respect to the atmosphere: The photosynthetic pathway used by plants is the primary determinant of the magnitude of the isotopic difference between atmospheric CO 2 , currently approximately −8.5‰ VPDB, and plant organic carbon. Plants using the C 3 photosynthetic pathway, which dominate most of the sites analyzed here, exhibit a wide FIORELLA ET AL.

10.1029/2020JG005862
3 of 20 range of isotopic compositions of −20‰ to ∼ −35‰ (Tipple & Pagani, 2007) that are quite distinct from atmospheric values. A common representation of carbon isotope discrimination at the leaf level relates two distinct kinetic fractionations with respect to atmospheric CO 2 , weighted by the concentration gradient between the atmosphere and that inside leaves: where a is the fractionation due to diffusion (4.4‰), and b is the fractionation associated with Rubisco carboxylation (27‰), c i is the intercellular (i.e., sub-stomatal) partial pressure of CO 2 and c a is atmospheric partial pressure of CO 2 (Farquhar et al., 1989). This expression is a simplified representation of carbon isotope fractionation and omits potential impacts of photorespiration (Farquhar et al., 1982;Seibt et al., 2008;Wingate et al., 2007) and smaller fractionations associated with molecular diffusion across the leaf boundary layer and dissolution of CO 2 in the leaf mesophyll (e.g., Farquhar et al., 1989;Suits et al., 2005).
In contrast to C 3 plants, carbon isotope fractionation in C 4 plants (many warm-weather herbaceous grasses and sedges) is both less variable and smaller in magnitude, with C 4 plant carbon isotope ratios typically ranging between −10 and −14‰. The primary difference between C 3 and C 4 fractionation arises from the first enzyme carboxylated during carbon uptake. Further details on carbon isotope fractionation during photosynthesis at the leaf-level and its response to environmental variability can be found in the literature (Cernusak et al., 2013;Farquhar et al., 1989;Tipple & Pagani, 2007).

Carbon Isotope Fluxes
Photosynthetic fractionation creates an isotopic signal that propagates through the ecosystem and can be used to probe ecosystem-scale processes. As described above, carbon assimilation selects for 12 C, which in isolation would leave the remaining CO 2 comparatively enriched in 13 C. Isotopic discrimination during photosynthesis is difficult to measure at ecosystem scales (e.g., Tissue et al., 2006;Wingate et al., 2010), and as a result, it is more common to estimate discrimination using a photosynthetic model for photosynthesis (Bowling et al., 2001;Knohl & Buchmann, 2005;Ogée et al., 2003;Wehr & Saleska, 2015). In turn, photosynthetic discrimination sets the baseline for post-photosynthetic plant metabolism and ecosystem respiration.
Yet, projecting leaf-scale controls on assimilatory and respiratory carbon isotope ratios to the ecosystem scale remains a challenge (e.g., Unger et al., 2010). Keeling plots, which rely on the covariance of δ 13 C-CO 2 and CO 2 mole fractions to estimate the isotope ratio of the net CO 2 flux, are the most common way to estimate the isotope ratio of ecosystem-scale respired CO 2 (Keeling, 1958;Pataki et al., 2003). Interpretations of carbon isotope ratios made from ecosystem-scale δ 13 C measurements can diverge from those that would be made with leaf-level measurements. For example, feedbacks between the land surface and the atmosphere can decouple ecosystem-level water use efficiency (e.g., assimilation/transpiration) from intrinsic water use efficiency recorded by leaf δ 13 C (Baldocchi & Bowling, 2003;Jarvis & McNaughton, 1986;Medlyn et al., 2017;Seibt et al., 2008). NEON observations of δ 13 C values in the atmosphere and in plant organic matter, in addition to new ecosystem modeling frameworks, provide the opportunity with harmonized data collection to test the applicability and generality of these relationships across scales.

Calibration Data
We implement a series of data screening and correction methods and use two different methods to calibrate observed ambient data to the VPDB scale. The first calibration method, modified from a method proposed by Bowling et al. (2003), estimates calibration parameters for each isotopologue ( 12 CO 2 and 13 CO 2 ) independently. The calibrated isotope ratios and CO 2 mole fractions are then calculated directly from the calibrated isotopologue mole fractions. The second method uses ordinary least squares (OLS) linear regression to generate an equation relating measured and known reference material isotope ratios to directly correct observed isotope ratios and CO 2 mole fractions without intermediate conversions to isotopologue mole fractions. While these methods are similar, two of the calibration gases often have similar δ 13 C values, and therefore, 12 CO 2 and 13 CO 2 mole fractions may be more broadly distributed through the measurement range than δ 13 C values, potentially yielding superior calibration from the Bowling et al. (2003) approach.
Internal instrument calibration parameters are corrected annually at the NEON Calibration, Validation, and Audit Laboratory in Boulder, which ensures that the average bias of the reported values for properly functioning instruments is reasonably small, but day-to-day instrumental drift can cause substantial deviations from these parameters. Calibration parameters are calculated using two calibration gas measurement periods bracketing a period (typically a day) of ambient atmosphere measurements. We impose a set of quality checks, common to both calibration approaches, on the calibration gas measurements that are distinct from the quality metrics calculated by eddy4R. First, calibration gas measurements are used only if known and measured CO 2 concentrations and δ 13 C values are all present. Second, calibration gas observations are used to calculate calibration parameters only if the measured CO 2 concentration of the calibration gas is within 10 ppm of the known CO 2 concentration and the variance of the CO 2 concentration over the 9 min measurement interval is less than 5 ppm. This filter helps to remove standard observations where the tank may have been empty or where there are other issues with the gas delivery system. Third, measured mean δ 13 C values are required to be within 3‰ of the known reference values. The choice of a 3‰ threshold is somewhat arbitrary and arises due to an observed gap in the distribution of the difference between measured and reference δ 13 C values between clearly erroneous and higher quality measurements. In addition, several sites were found to have incorrect reference gas values in the currently provisioned files; we have provided corrected values in the Table S1. This filter generally removes only a very small fraction of reference gas measurements, which most typically occur in the first day or two of measurements after an instrument has been inactive for an extended period. The code for correcting and calibrating the NEON carbon isotope data described here is available as an R package, NEONiso, that can be installed from GitHub (https://doi.org/10.5281/zenodo.3836875). Version 0.3 of the NEONiso package was used to produce the calibrated data presented here. The NEON data comprising this analysis are very similar to those contained within the 2021 EC data release issued by NEON (RELEASE-2021; 2021); deviations from this data release are described in the Tables S1 and S2.

Gain-And-Offset Method Calibrating Individual Isotopologues
The Bowling et al. (2003) approach seeks to correct 12 CO 2 and 13 CO 2 mole fractions independently and uses these values to calculate calibrated CO 2 and δ 13 C values. Mole fractions for measured and reference 12 CO 2 and 13 CO 2 are calculated from measured and reference CO 2 and δ 13 C values using the definition of delta notation and mass balance considerations: Equation 1 can be rearranged for the ratio of 13 CO 2 and 12 CO 2 in the measured gases as: Considering the total CO 2 mole fraction as a sum of all possible isotopologues yields: where f (0.00474) is the fraction of CO 2 comprised of isotopologues other than 12 C 16 O 2 or 13 C 16 O 2 (e.g., 14 CO 2 , 12 C 16 O 18 O). In theory, f varies with the isotope ratio of the measured gas but these variations are negligibly small in practice and therefore f is treated as constant (Griffis et al., 2004). Equation 4 and 5 can be rearranged to give expressions for [ 12 CO 2 ] and [ 13 CO 2 ]: Gain and offset factors to correct measured 12 CO 2 and 13 CO 2 mole fractions based on reference values are calculated using measurements of two different reference gases. For reference gases A and B, gain factors (G i ) and offsets (O i ) are calculated as: where X B,ref and X A,ref are the known mole fractions for the reference gases and X B,meas and X A,meas are the measured mole fractions of isotopologue i (either 12 CO 2 or 13 CO 2 ), respectively.
Calibrated mole fractions of 12 CO 2 and 13 CO 2 are then calculated from these gain and offset functions as: While the initial implementation in Bowling et al. (2003) used only two calibration gases, using the third as a check on the calibration, we have estimated the gain and offset parameters using linear regression of measured and reference values of 12 CO 2 and 13 CO 2 using all three reference gases from two consecutive reference gas measurement periods. If the R 2 value of the regression is ≥ 0.95 for both isotopologues, the period of ambient measurements bracketed by the data used in this regression are calculated based on these gain and offset parameters. Using the calibrated 12 CO 2 and 13 CO 2 mole fractions from Equation 10, δ 13 C values are calculated with Equation 1, and calibrated CO 2 values are calculated by rearranging Equation 5: Where the R 2 of these regressions was below 0.95, ambient measurements were excluded from further analysis.

Direct Calibration of δ 13 C Values and CO 2 Mole Fractions Using Linear Regression
We have also developed a calibration routine using linear regression to directly calibrate measured isotope ratios to the reference VPDB scale. Calibration gas measurements from two adjacent periods are used to estimate the slopes and intercepts of lines to correct δ 13 C and CO 2 values to the known values. If all candidate data points passed the initial quality checks, six data points are used for each linear regression (two for each reference material, spread one day apart). If the R 2 value of the regression is ≥ 0.95, this regression is used to calibrate the period of ambient measurements bracketed by the data used in this regression using the estimated slope (m) and intercept (b) values: where the R 2 of either regression was below 0.95, ambient measurements were excluded from further analysis.

Post-Calibration Data Processing
Despite the QA/QC filtering steps above, a small fraction (<0.01% for both calibration methods) of calibrated δ 13 C values were outside of a plausible range, which we defined as −25‰ to 0‰ and were removed from subsequent analysis. Impulse spikes in the δ 13 C and CO 2 measurements were removed using a median-absolute-deviation filter (Brock, 1986;Starkenburg et al., 2016). Additionally, a few sites exhibited an FIORELLA ET AL.

10.1029/2020JG005862
6 of 20 apparent step change in isotope ratios that persisted for at least a month and could not be tied to an obvious calibration issue. The cause of these issues was not immediately apparent but could be related to instrument or sampling manifold performance issues. Visual inspection of calibrated timeseries led us to remove 8 months of data from further analysis: (a) May to September 2019 for UNDE, (b) July 2019 for SRER, and (c) August and September 2018 for TEAK. After these months were removed, the processed dataset consisted of 1,654 months from 47 NEON sites.
Calibration uncertainties are reported as both differences (to reflect bias) and root-mean-square errors (RMSE, to reflect uncertainty) using the predicted calibrated δ 13 C and CO 2 values relative to the known reference δ 13 C and CO 2 values. To generate these error estimates, a leave-one-out analysis was used, with calibration routines applied to two of the three standards and leaving the third to estimate the error. This procedure was repeated three times, leaving each standard out in turn, and is analogous to a jackknife style error estimation procedure. Errors for all three permutations were pooled together to generate a distribution of error estimates for each site. This procedure likely overestimates the true calibration error, as it only uses information from two of the reference gases while our final calibrated data used all three reference gases, and thus were constrained by a greater amount of information.
For a standardized estimate of the difference between calibrated and uncalibrated values across sites, we also calculated a root-mean-squared difference for each ambient measurement point. Seasonal comparisons are based on meteorological seasons: summer is defined as June, July and August while winter is defined as December, January, and February. For analyses on diurnal timescales, we approximate night as 10 p.m. to 5 a.m. local time and day as 10 a.m. to 5 p.m. local time. Site elevations and mean annual precipitation and temperature are derived from NEON-provided metadata (Table 1), while analyses using an aridity index, approximated by precipitation over potential evapotranspiration (P/PET) are extracted from the GCIAR Consortium for Spatial Information Global Aridity Index dataset (version 2), using the pixel value corresponding to the NEON tower location (Zomer et al., 2008). NEON sites were divided into forest, grassland, shrubland, tundra, or agricultural classes based on the dominant NLCD land class reported by NEON for their tower footprint plots (Table 1). We also estimate the δ 13 C of respired CO 2 (hereafter δ 13 C R ) to gain insight into large-scale controls on ecosystem-level carbon isotope fractionation. To estimate δ 13 C R we use a Keeling style plot by regressing the product of the observed CO 2 mole fraction and δ 13 C value against the measured CO 2 mole fraction in the nighttime data during JJA (Miller & Tans, 2003). Site-mean summer δ 13 C R values are compared against mean summer maximum daily vapor pressure deficit (VPD) to explore the ecosystem response to evaporative demand across ecosystems. Sites were classified as either primarily C 3 or mixed C 3 /C 4 ecosystems based on whether C 4 species were present within the flux-tower survey plots in NEON's "Plant presence and percent cover" data sets (data product DP1.10058.001; NEON., 2020).

Performance of Calibration Methods
Calibrated isotope ratios and CO 2 concentrations across the NEON network are similar for both calibration methods (Figure 1). For δ 13 C, the between-calibration differences were <1‰ for >99.9% of data points and <0.25‰ for 98% of points. Deviations in δ 13 C values away from the 1:1 line suggest rare periods where one calibration method deviates strongly from the other. Differences in CO 2 mole fraction are also generally small across the two calibration methods, with both calibrations yielding CO 2 mole fractions that were always within 1 ppm of each other (Figure 1b).
Despite similar results, the gain-and-offset approach produced lower-error calibrations than directly calibrating δ 13 C values and CO 2 mole fractions. For 31 of 47 NEON sites reporting atmospheric carbon isotope data, the gain-and-offset approach produced lower δ 13 C error estimates than direct calibration of δ 13 C values (Figures 2a and 2c), but with lower calibration quality for CO 2 mole fractions (Figures 2b and 2d). In addition, high CO 2 RMSE values were observed for a few sites where there were issues with reference gas delivery, most notably SCBI (Figure 2d). Across all sites, the median bias for δ 13 C (CO 2 ) was 0.05‰ (0.04 ppm) and the uncertainty, quantified by RMSE, for δ 13 C (CO 2 ) was 0.40‰ (0.26 ppm) for the gain-and-offset calibration. Direct calibration of δ 13 C and CO 2 values yielded median bias estimates of 0.11‰ (0.01 ppm) for δ 13 C (CO 2 ), and RMSE values for δ 13 C and CO 2 of 0.42‰ (0.15 ppm). For the gain-and-offset calibration, FIORELLA ET AL. 72% and 15% of sites had mean RMSE better than 0.5‰ and 0.25‰ for δ 13 C values respectively, whereas these thresholds were met for 66% and 21% of sites for the direct δ 13 C and CO 2 calibration approach. Larger uncertainties in the calibration of δ 13 C, but not CO 2 , are observed at most sites for the direct δ 13 C and CO 2 calibration approach compared to the gain-and-offset calibration approach due to an increase in high-δ 13 Cerror calibrations ( Figure S1), driving a larger range of RMSE values ( Figure S2). The most likely explanation for this pattern is that for many sites, two of the reference gases have similar δ 13 C values but different CO 2 mole fractions; in this setup, estimation of the third reference gas δ 13 C value would have a high error due to a large uncertainty in the calibration slope. Estimation of total CO 2 , or of 12 CO 2 and 13 CO 2 in the gain-and-offset method, is not subject to the same limitation as there is a larger variation in these quantities despite the similarity of the δ 13 C values. As these calibration methods yielded similar performance, we focus our results on calibrated carbon isotope data using the gain-and-offset calibration.
The magnitude of the correction applied in calibrating necessary to yield calibrated δ 13 C and CO 2 mole fractions was highly site specific, underscoring the need to apply a consistent calibration approach to each site of a multisite network ( Figure S3). Median root-mean-square differences for δ 13 C varied from 1.07‰ (UNDE) to 0.14‰ (YELL), with occasional larger excursions of up to 2‰ observed for some sites ( Figure S3a). The range of median root-mean-square differences for CO 2 was 2.9 ppm (SCBI) to <0.1 ppm (ORNL), with occasional deviations greater than 5 ppm ( Figure S3b). These differences are likely a reflection of varying drift characteristics of each instrument across time, and apparently poor internal calibration of CO 2 mole fractions or mismatched calibration metadata for a handful of NEON sites (Table S1). The magnitudes of FIORELLA ET AL.
10.1029/2020JG005862 9 of 20 these changes with calibration and their variations across time and space demonstrate that necessity of applying these repeated recalibrations. Fortunately, correcting both δ 13 C values and CO 2 mole fractions using a consistent calibration framework is possible, which would enable robust cross-site comparisons.

Cross-Site Annual Patterns
The post-calibration data reveal cross-site responses to environmental drivers. Mean δ 13 C isotope ratios across all available data for each site exhibit negative correlations with mean annual precipitation (Figure 3a), aridity ( Figure 3b) and mean annual temperature (Figure 3c), but a positive correlation with site elevation (Figure 3d). Isotopic discrimination during assimilation is known to be higher when plants are not water stressed and stomata are open, promoting higher c i /c a ratios and lower 13 C values (e.g., Cernusak et al., 2013), which likely underlies the observed relationships with mean annual precipitation and aridity. Reduced carbon isotope discrimination has also been observed at high elevations in plant samples (e.g., Körner et al., 1988), consistent with our observations of higher mean δ 13 C in atmospheric CO 2 at high elevations ( Figure 3d). The mechanisms responsible for decreased carbon discrimination at high elevation remain debated; for example, both increases in plant photosynthetic capacity and changes in water vapor and oxygen partial pressures perhaps contribute (e.g., Wang et al., 2017). Data shown in Figure 4 are from the lowest tower level at each site, but qualitatively equivalent patterns are observed when considering the highest tower level at each site.
Differences in δ 13 C ( Figure 4) and CO 2 ( Figure 5) mole fractions across the tower depth were most pronounced in forests, less pronounced at agricultural sites, and nearly absent in shrubland, grassland, and tundra sites. Observed differences in annual median δ 13 C between top and bottom measurement levels ranged from −0.06‰ (OSBS) to 0.94‰ (TEAK) at forest sites, and from 0.02‰ (KONA) to 0.34‰ (LAJA) at agricultural sites. In contrast, only one tundra site (NIWO), one grassland site (CPER), and no shrubland sites had a difference in median δ 13 C larger in magnitude than 0.15‰. Similarly, vertical differences in CO 2 mole fraction were also greater in forest and agricultural sites than in the shrubland, grassland, or tundra sites. Top-to-bottom differences in annual median CO 2 values ranged from −21.6 ppm (TEAK) to 0.6 ppm (OSBS) at the forest sites, and −0.2 ppm (KONA) to −20.1 ppm (LAJA) at the agricultural sites. The larger gradients at forest sites are likely attributable to both higher carbon uptake (e.g., JERC) and de-coupling of air exchange between canopy levels (e.g., WREF), whereas the gradients at agricultural sites are likely only attributable to higher carbon uptake. Across grassland, shrubland, and tundra sites, only BARR (3.4 ppm), NIWO (−1.7 ppm) and CPER (−1.2 ppm) exhibited annual CO 2 differences from tower top to bottom larger than 1 ppm. Further investigation into 13 C variations across sites, and how they may be driven by assimilation and respiration rates (section 2), can be further evaluated as a more extensive and continuous CO 2 EC dataset becomes available.

Temporal Variation in Vertical Profiles
Vertical differences observed in δ 13 C values and CO 2 mole fractions vary across seasons and time of day.
Median top-to-bottom profile differences are larger in magnitude in summer than winter for 30 of 47 NEON sites for δ 13 C ( Figure S4) and 35 of 47 sites for CO 2 ( Figure S5). As was observed at the annual scale, the largest magnitude differences through the tower depth are observed during summer at forest sites. The average δ 13 C top-to-bottom difference in summer was 0.70‰ for forested sites, but 0.02‰, <0.01‰, −0.05‰, 0.02‰ for agricultural, shrubland, tundra, and grassland sites, respectively. δ 13 C differences between tower top and bottom were smaller in magnitude during winter at all site classes except for agricultural sites: FIORELLA ET AL. forests, shrublands, tundra, and grassland sites exhibited winter differences of 0.23, 0.03, 0.01, and −0.05‰ respectively, while agricultural sites exhibited a winter δ 13 C difference of 0.22‰.
Similar patterns are observed for CO 2 mole fractions ( Figure S5). Average CO 2 concentrations at the top of the tower in forested sites are 18.8 ppm lower than CO 2 concentrations at the lowest tower level during summer, but this difference is only 2.6 ppm and 1.1 ppm for agricultural and grassland sites, respectively. Summertime CO 2 concentration gradients are reversed but lower magnitude at shrubland and tundra sites, where CO 2 concentrations are higher at the tower top by 0.2 and 1.5 ppm, respectively, than at the tower base. Winter CO 2 differences between top and bottom tower intakes are lower in magnitude for forest (−4.6 ppm), shrubland (−0.3 ppm), tundra (−0.3 ppm), and grassland (0.5 ppm) sites, but higher than summer values for agricultural sites (−9.2 ppm).
The larger magnitudes of CO 2 and δ 13 C gradients in summer relative to winter likely reflect that photosynthetic and respiratory fluxes are generally highest during the growing season, which is meteorological summer for most of these sites. Multiple factors may contribute to the larger differences observed in forests, including: (1) photosynthetic and respiratory fluxes are high in forests (Beer et al., 2010;Bonan, 2008;Buchmann & Schultze, 1999), and (2) forest canopies may be partially decoupled from the overlying atmosphere, which inhibits mixing that would reduce the gradient (Aron et al., 2019;Jarvis & McNaughton, 1986;Rastogi et al., 2018). In contrast, winter patterns for the agricultural sites are primarily driven by a single site in this category, LAJA, which appears to have a dominant growing season in meteorological winter.
Substantial summer variations in δ 13 C ( Figure 6) values and CO 2 (Figure 7) mole fractions differences are also observed on diurnal timescales. Most sites exhibit marked decreases in near-surface δ 13 C and increases in CO 2 mole fractions at night, and a more uniform profile during the day. All sites but five (BARR, TOOL, KONZ, OAES, and MLBS) exhibit δ 13 C gradients that are larger in magnitude at night than during day FIORELLA ET AL.

10.1029/2020JG005862
13 of 20 ( Figure 6), and only BARR, TOOL, and OAES exhibit a larger magnitude CO 2 gradient during the day compared to during night. All ecosystem types exhibit larger δ 13 C top-to-bottom differences at night, on average, though this pattern is most pronounced in forest ecosystems (1.37‰ at night, 0.29‰ during the day). Differences in CO 2 mole fraction (Figure 7) are most pronounced in forest (−39.9 ppm from top to bottom at night, −7.1 ppm during the day) and agricultural ecosystems (−41.4 ppm at night, 2.1 ppm during the day). Larger nighttime gradients are the product of reduced atmospheric mixing at night, which allows the buildup of respired CO 2 near the surface and allows estimation of the δ 13 C value of the respiratory flux using Keeling plots (e.g., Pataki et al., 2003). The comparatively small change in δ 13 C gradients, but large change in CO 2 mole fraction gradients observed at agricultural sites likely results from the substantial presence of C 4 vegetation at all four agricultural sites (LAJA, STER, KONA, and DSNY large increases in respiratory CO 2 , as they are highly productive, but small decreases in δ 13 C associated with the respiratory flux since the isotopic composition of organic carbon produced by C 4 vegetation is closer to the atmospheric value . Strong variations in the estimated δ 13 C value of the summer respiration flux were observed across the NEON network ( Figure 8), with estimated JJA δ 13 C R values ranging from −29.4‰ (NIWO) to −16.6‰ (LAJA). Sites with primarily C 3 vegetation had lower mean δ 13 C R values (−26.9 ± 0.9‰) than sites with mixed C 3 /C 4 vegetation (−22.3 ± 2.9‰), with the larger standard deviation in mixed C 3 /C 4 sites likely reflecting variations in the C 4 abundance and role in the local carbon cycle across these sites. JJA δ 13 C R values at C 3 -dominated sites also exhibited a strong dependence on VPD, with higher atmospheric demand corresponding to a FIORELLA ET AL. higher value of δ 13 C R (Figure 8, slope = 0.58 ± 0.18‰/kPa, r = 0.516, p = 0.002). In contrast, no relationship between δ 13 C R values at mixed C 3 /C 4 sites and VPD was observed.

Conclusions
We have outlined two potential calibration strategies for CO 2 isotope ratio measurements made as part of NEON EC storage exchange measurement system. We find that independent calibration of 12 CO 2 and 13 CO 2 (after Bowling et al., 2003) using linear regression produces a lower error calibration than a direct calibration of δ 13 C and CO 2 values, but at a slight cost to CO 2 mixing ratio precision. Wen  similar conclusion with an earlier generation of the Picarro carbon isotope analyzer. Furthermore, we have demonstrated the importance of applying a consistent calibration technique across all NEON sites, as: (a) deviations between uncalibrated and the calibrated, drift-corrected δ 13 C and CO 2 values were strongly site dependent, and (b) the removed errors were often comparable in magnitude to the signal of interest (Figures 2, S1-S3). In addition to the functionality demonstrated here, the NEONiso R package also contains functions to produce diagnostic plots of calibration data and ambient δ 13 C-CO 2 and CO 2 mole fractions. Future developments to the NEONiso R package, as well as anticipated operational improvements at NEON field sites with time, will likely further improve the quality of the calibration, as well as expand the range of integrative analyses that can be performed with NEON data products.
We have documented broad network-wide patterns in this work, demonstrating the ability of the post-calibration NEON CO 2 isotope ratios to capture meaningful macroscale biogeochemical and ecological patterns.
Mean δ 13 C values at NEON sites appear to vary with mean annual precipitation, aridity, mean annual temperature, and site elevation (Figure 3). These findings echo earlier works based principally on isotopic analysis of leaf material (Diefendorf et al., 2010;Farquhar et al., 1989;Kohn, 2010;Körner et al., 1988Körner et al., , 1991. Gradients within the tower observations vary strongly with ecosystem type, season, and time of day (Figures 4-7). Gradients are largest within forested ecosystems, particularly during the summer, when NEE is higher and more variable, and at night, when the overlying atmosphere is stable. The calibrated data products described here will enable much more advanced analyses than have been previously possible by merging these data with the larger catalog of ecological and meteorological data collected by NEON. For example, while some studies have shown a strong response of the δ 13 C value of the respiration flux to VPD (Bowling et al., 2002), other studies have shown that this relationship is highly site dependent (e.g., Raczka et al., 2017). Standardized data collection protocols over a large range of ecosystem types and climates may help clarify this relationship, and the methods proposed here could be used for standardization of future continental-scale observations of atmospheric stable isotope ratio measurements. Similarly, canopy conductance is often estimated using a re-arrangement of the Penman-Monteith equation and latent heat fluxes (e.g., Novick et al., 2016), but the addition of carbon isotope ratios allows estimation of canopy conductance from the estimated assimilatory flux and δ 13 C value (Griffis et al., 2004;Knohl & Buchmann, 2005;Wehr & Saleska, 2015). As NEON data catalogs of EC and atmospheric isotope ratio data continue to grow, it will be possible to constrain isotope ratios of individual carbon fluxes (e.g., Bowling et al., 2001Bowling et al., , 2014 and improve our understanding of individual processes driving NEE and their representation in ecosystem and land surface models.

Data Availability Statement
Code to calibrate NEON carbon isotope data are available as the NEONiso R package from GitHub: https:// doi.org/10.5281/zenodo.3836875. The CGIAR Global Aridity Dataset are available for download from: https://figshare.com/articles/Global_Aridity_Index_and_Potential_Evapotranspiration_ET0_Climate_Da-tabase_v2/7504448/3. NEON data are available from the NEON data portal: https://data.neonscience.org. The National Ecological Observatory Network program is sponsored by the National Science Foundation and operated under cooperative agreement by Battelle Memorial Institute.
FIORELLA ET AL.

10.1029/2020JG005862
17 of 20 Figure 8. Variation in the estimated JJA δ 13 C value of ecosystem respiration (δ 13 C R ) with mean daily maximum VPD. A positive correlation between δ 13 C R and VPD was observed in sites dominated by plants using the C 3 photosynthetic pathway. The best fit line between δ 13 C R and VPD at C 3 sites is shown as a black line, with gray shading indicating the 95% confidence interval of the regression. In contrast, sites that also had plants using the C 4 photosynthetic pathway tended to have higher δ 13 C R that were not correlated with VPD.