The Berkeley Environmental Air-quality and CO₂ Network: field calibrations of sensor temperature dependence and assessment of network scale CO₂ accuracy

The majority of global anthropogenic CO₂ emissions originate in cities. We have proposed that dense networks are a strategy for tracking changes to the processes contributing to urban CO₂ emissions and suggested that a network with ∼ 2 km measurement spacing and ∼ 1 ppm node-to-node precision would be effective at constraining point, line, and area sources within cities. Here, we report on an assessment of the accuracy of the Berkeley Environmental Air-quality and CO₂ Network (BEACO₂N) CO₂ measurements over several years of deployment. We describe a new procedure for improving network accuracy that accounts for and corrects the temperature-dependent zero offset of the Vaisala CarboCap GMP343 CO₂ sensors used. With this correction we show that a total error of 1.6 ppm or less can be achieved for networks that have a calibrated reference location and 3.6 ppm for networks without a calibrated reference.

1 Introduction

The atmosphere has warmed approximately 1 ± 0.2 ^∘C since pre-industrial times, which is unequivocally due to anthropogenic emissions of CO₂ and other greenhouse gases (GHGs) (IPCC, 2021). Global initiatives are needed to limit warming to 1.5 ^∘C by achieving net zero GHG emissions by 2050 and a 45 % emissions decline from 2010 levels by 2030 (Rogelj et al., 2021). As over 70 % of global anthropogenic CO₂ emissions originate from cities (United Nations, 2011), effective CO₂ monitoring strategies in urban regions are needed to assess progress toward emissions commitments.

Monitoring trends in CO₂ emissions by tracking ambient CO₂ in urban environments is challenging because of the large diversity of emissions sources, complex spatial and temporal patterns of emission rates, varied topography, and the effects of meteorology on the observed concentrations (e.g., Vardoulakis et al., 2003; Lateb et al., 2016). As a result, most cities rely exclusively on economics and social data and do not check whether their reported emissions match the observed CO₂ enhancements in the air over their city. To date, most efforts to assess CO₂ emissions from cities have relied upon a small number of high-cost CO₂ instruments that provide precise and accurate representations of regional signals. Other approaches include use of correlations between CO₂ and other gases, measurements of ¹⁴C in annual grasses, and use of satellite column CO₂ observations such as from OCO-2 (e.g., Pataki et al., 2003, 2006; Riley et al., 2008; Thompson et al., 2009; Kort et al., 2013; Andrews et al., 2014; Fu et al., 2019; Ye et al., 2020). Most of these efforts have used as a target metric an annual average of fossil fuel-related CO₂ emissions from an entire city (e.g., McKain et al., 2012; Kort et al., 2013; Bréon et al., 2015; Verhulst et al., 2017). Simultaneous measurements of CO and ¹⁴CO₂ have also provided information about sector-specific emission sources (Turnbull et al., 2015). Other methods of evaluating urban emissions have relied on emissions inventories (e.g., Gurney et al., 2009; Gately et al., 2013, 2017). These emissions inventories are frequently applied to inverse modeling approaches in combination with either short-term mobile measurements or a small number of long-term measurement sites to extract regional emissions (e.g., Brondfield et al., 2012; Sargent et al., 2018; Nathan et al., 2018; Turnbull et al., 2019). Several studies have also combined a network of CO₂ observations with inverse modeling approaches to evaluate the accuracy of emissions inventories and CO₂ sources (e.g., Lauvaux et al., 2016, 2020).

We are pursuing a distinct approach aimed at process-level understanding of the components of an urban emissions inventory. To do so, we are developing tools for the deployment of spatially dense networks of CO₂ measurements, in combination with gases and aerosols that are co-emitted and that affect air quality. The result is an ability to map emissions with ∼ 1 km or “neighborhood scale” fidelity. The Berkeley Environmental Air-quality and CO₂ Network (BEACO₂N) (Turner et al., 2020; Kim et al., 2018; Shusterman et al., 2018, 2016; Turner et al., 2016) is our platform for research and development of tools for dense networks. New deployments in Glasgow, Scotland and Los Angeles, California are bringing new collaborators and experience in different cities to the project. BEACO₂N has been operating since 2012 in the San Francisco Bay area and consists of over 70 nodes separated by approximately 2 km (Fig. 1). The nodes incorporate commercially available, low-cost sensors for measuring CO, NO, NO₂, O₃, particulates, and CO₂.

Figure 1 Map of all Bay Area BEACO₂N sites (small red dots), BEACO₂N sites discussed in this work (large blue dots), and the Richmond Field Station (star).

Turner et al. (2016) assessed the performance of a hypothetical BEACO₂N-like observing system coupled to an inverse model and demonstrated that a random measurement uncertainty of 1 ppm between nodes was adequate to meaningfully constrain CO₂ emissions from a point, line, or area source of 147, 45, and 9 t C h⁻¹, respectively. With a 1 ppm mismatch error, weekly CO₂ emissions in the San Francisco Bay area could be estimated to within 5 % error. In this paper we describe advances in our approach to maintaining stable, multiyear comparability among BEACO₂N nodes in a city and evaluate the accuracy achieved with these new procedures. Our emphasis in the revised approach to sensor accuracy is on tracking and correcting the temperature dependence of the Vaisala CarboCap GMP343 CO₂ instruments. We present the development and evaluation of the methods using observations from the San Francisco Bay Area BEACO₂N deployment and then apply these ideas to the BEACO₂N network in Houston, Texas.

2 Development of a CO₂ field calibration method for Vaisala temperature dependence

The efficacy of a network of a large number of low-cost nondispersive infrared (NDIR) CO₂ sensors to evaluate CO₂ emissions has been previously discussed (Shusterman et al., 2016; Turner et al., 2016; Martin et al., 2017; Shusterman et al., 2018; Müller et al., 2020). Martin et al. (2017) showed that after correcting six SenseAir K30 carbon dioxide NDIR sensors (with off-the-shelf reported errors of 5–20 ppm) for environmental variables, the median root mean square error could be reduced to below 2 ppm, making the sensors potentially useful for ambient air-quality monitoring. Recently, Müller et al. (2020) evaluated the potential applications of a low-cost CO₂ NDIR sensor network for resolving site-specific CO₂ signals in Switzerland. The calibration method of Müller et al. (2020) involved laboratory chamber calibrations of over 300 low-cost NDIR CO₂ sensors and ambient colocation with a reference instrument prior to deployment, as well as regular monitoring and drift correction during a 2-year deployment period. Shusterman et al. (2016) developed an in situ method for calibrating and correcting for individual instrument biases and temporal drifts of the Vaisala CarboCap GMP343 CO₂ instruments deployed in the BEACO₂N nodes. Using this method, Shusterman et al. (2018) demonstrated that the BEACO₂N network could provide highly sensitive detection of changes to traffic emissions at a scale relevant to policy concerns. Shusterman et al. (2018) also illustrated the efficacy of the BEACO₂N network in showing both regional CO₂ emissions and local CO₂ enhancements at the scale of a single neighborhood. In an analysis of the BEACO₂N observations for 6 weeks before and after the COVID-19 shutdown, Turner et al. (2020) showed that a 25 % change in emissions is easily derived by an inverse model and that hourly variations in emissions can be inferred.

The use of a large number of low-cost CO₂ sensors introduces challenges regarding accuracy and inconsistent behavior between instruments that often requires labor-intensive regular calibration, data correction and filtering, and validation with comparison to a smaller number of frequently calibrated high-accuracy instruments. In particular, the low-cost NDIR absorption sensor used in each BEACO₂N node (Vaisala CarboCap GMP343) is susceptible to temporal drift and fluctuations due to environmental variables that present challenges to achieving a goal of 1 ppm network error (van Leeuwen, 2010; Shusterman et al., 2016). Correction of the Vaisala CarboCap GMP343 instruments (Vaisala, hereafter) for changes in pressure, temperature, and humidity is required for accurate measurements (Vaisala, 2013). The typical correction for pressure and temperature accounts for changes in the number density of CO₂ according to the ideal gas law (van Leeuwen, 2010; Vaisala, 2013; Shusterman et al., 2016). The humidity effect on measured CO₂ is accounted for by considering the dilution effect of water vapor according to Dalton's law of partial pressures (van Leeuwen, 2010; Vaisala, 2013; Shusterman et al., 2016). However, even after accounting for these factors, reported corrected CO₂ concentrations for the Vaisala instrument have been observed to exhibit a strong temperature dependence of up to 1 ppm ^∘C⁻¹ (van Leeuwen, 2010). Using a laboratory calibration procedure, van Leeuwen (2010) found that a linear correction was necessary to account for the residual temperature dependence. However, correcting for the temperature dependence using lab calibrations is labor intensive, as the temperature dependence is unique for each Vaisala sensor. Regular laboratory temperature calibration would also be required to account for temporal variations in the temperature correction as sensors age. For a high-density urban network like BEACO₂N, this would require substantial time investment by trained personnel. The associated high labor costs defeat the purpose of using low-cost sensors. In situ field calibration of the Vaisala sensors thus presents a more attractive method for correcting for the temperature dependence of the CO₂ measurements.

2.1 BEACO₂N network

The Berkeley Environmental Air-quality and CO₂ Network (BEACO₂N) Bay Area deployment currently consists of 73 nodes spaced at ∼ 2 km intervals with locations in Alameda, San Francisco, Contra Costa, Sonoma, Sacramento, and Solano counties. A full description of a BEACO₂N node can be found in Kim et al. (2018). Briefly, each node contains a nondispersive infrared Vaisala CarboCap GMP343 CO₂ sensor, along with a Shinyei PPD42NS nephelometric particulate matter sensor and several Alphasense electrochemical sensors for measuring CO, NO, NO₂, and O₃ (CO-B4, NO-B4, either NO2-B42F or NO2-B43F, and either Ox-B421 or Ox-B431). The most recent version adds a Plantower PMS 5003 aerosol sensor. Sensors are assembled into compact, weatherproof enclosures with air flow through the enclosure provided by two 30 mm fans. Data are compiled with a Raspberry Pi microprocessor and an Adafruit Metro Mini microcontroller. Data are acquired every 5 or 10 s and are transferred to a central server via an Ethernet or Wi-Fi connection. Observations are posted on the BEACO₂N website within a few hours of measurement time (http://beacon.berkeley.edu, last access: 14 June 2021).

The Vaisala CarboCap GMP343 instrument uses pulsed light from a filament lamp, which is reflected and refocused on an IR detector located behind a Fabry-Perot Interferometer (FPI). The FPI is electrically tuned so that its passband corresponds to either the absorption wavelength of CO₂ or a reference band (Vaisala, 2013). The calibration procedure for the Vaisala CarboCap GMP343 CO₂ sensor is as outlined in Shusterman et al. (2016, 2018). Briefly, deployed Vaisala sensors operate with the internal relative humidity (RH), temperature, and pressure compensation set to “off” and the oxygen correction set to “on”, with oxygen input as 20.95 %. A post hoc multiplicative scale factor is applied to convert the raw CO₂ outputs to the mole fraction of CO₂ that would be measured if the observed air parcel were dried and brought to standard temperature and pressure (STP) ([CO₂]_STP). Raw CO₂ data are adjusted using temperature (T) measured by the internal thermometer of the Vaisala. Water vapor pressure ($P_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$) and air pressure (P_tot) are obtained from the pressure and dew point temperature measured inside each node enclosure by a Bosch Sensortec Adafruit BME280 sensor. The [CO₂]_STP is then adjusted to account for temporal drift in the instrument “zero” by comparing the background signal of the Vaisala CO₂ measurement at each node to a reference Picarro G2301 system, located at the Richmond Field Station in Richmond, CA (Fig. 1). A moving 3-week window of the 10th percentile of Vaisala CarboCap CO₂ data (^Vaisala[CO₂]_10 %) is generated and compared with the 10th percentile of the reference Picarro instrument (^Picarro[CO₂]_10 %). The difference between ^Vaisala[CO₂]_10 % and ^Picarro[CO₂]_10 % is used to define the offset of the Vaisala instrument $\left([{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{offset}}^{\mathrm{drift}}\right)$. A linear correlation between $\left([{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{offset}}^{\mathrm{drift}}\right)$ and time is generated and used to calculate the drift-corrected CO₂ data, $\left([{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{corrected}}^{\mathrm{drift}}\right)$ (Eqs. 1–2):

$$\begin{array}{}\text{(1)}& {\displaystyle }\left([{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{offset}}^{\mathrm{drift}}\right)=m_t\times \text{days}+b\text{(2)}& {\displaystyle }[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{corrected}}^{\mathrm{drift}}=[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{STP}}-m_t\times \text{days}-b,\end{array}$$

where m_t is the temporal drift (ppm d⁻¹) and b is a constant atemporal bias.

2.2 Picarro reference instrument

A “supersite” with reference grade instruments is operated within the BEACO₂N Bay Area network to provide a reference for the network calibration. Instruments are installed within a temperature-controlled instrument shelter at the U.C. Berkeley Richmond Field Station. Measurements include basic meteorology, NO_x (Thermo 42CTL with a molybdenum NO₂ to NO converter), O₃ (Teledyne/API T400), CO₂, CH₄, and CO (Picarro G2401 cavity ring down analyzer). Air is sampled through Teflon tubes mounted to a small tower affixed to the trailer roof, for a combined height of 6 meters above the ground. The colocated BEACO₂N node is attached outside of the trailer to the same tower.

The NO_x and Picarro analyzer calibrations are checked against reference gases every 2 to 3 weeks. The reference gas cylinders for NO_x, CO, and CH₄ are Certified Standard grade from Praxair, and for CO₂ are from the NOAA Global Monitoring Laboratory (two levels: 403.61 and 687.47 ppm). The Picarro checks are made by flowing the sequence of references gases into a tee at the inlet of the instrument for 15 min per step. The sequence of steps is performed twice during a check. The flow rate is set to be larger than the instrument sample flow (0.4 L min⁻¹) to overflow the inlet. The O₃ analyzer is checked against a photometric calibrator (Teledyne/API 703E).

2.3 Identification of a temperature-dependent error in Vaisala measurements

There exists an additional temperature dependence among the Vaisala CarboCap GMP343 instruments that varies between instruments. The temperature dependence was first identified from observations of CO₂ diurnal cycles at certain Bay Area BEACO₂N sites that were out of phase or larger in magnitude than the diurnal cycles at nearby nodes or measured by the Picarro. The presence of a temperature dependence in suspect Vaisala instruments was confirmed by examining the relationship between temperature in the node and the difference between baseline CO₂ signals measured by the Vaisala and the Picarro reference instrument.

Diurnal cycles of urban CO₂ typically exhibit a daily maximum at night or mid-morning (depending on influence from traffic emissions) due to mixing in a shallow nighttime planetary boundary layer (PBL), and reach a minimum during the day as PBL height increases and vegetation takes up CO₂ (Idso et al., 2002; Coutts et al., 2007; Turnbull et al., 2015; Shusterman et al., 2016). The presence of an additional temperature dependence in the Vaisala CO₂ instrument is particularly pronounced and obvious in the measurements obtained with the sensor located at the East Bay Municipal Utility District (EBMUD) BEACO₂N site during 2020 (Fig. 2). The magnitude of the diurnal cycle at EBMUD is larger and out of phase with the Picarro reference instrument (Fig. 2a). The result of this temperature dependence at EBMUD (Fig. 2c) is a diurnal cycle that peaks midday (Fig. 2b). Figure 2b compares the CO₂ diurnal cycle at EBMUD with the nearby urban site Laney College (Fig. 1). In contrast to EBMUD, Laney College exhibits a daily maximum at mid-morning–a pattern more consistent with typical urban CO₂ behavior (Idso et al., 2002; Coutts et al., 2007; Turnbull et al., 2015; Shusterman et al., 2016).

Figure 2 (a) CO₂ mixing ratios from April 2020 at EBMUD and measured with a Picarro instrument at the Richmond Field Station supersite. (b) EBMUD 2020 diurnal cycle compared with Laney College. (c) Temperature dependence of the CO₂ signal at EBMUD.

The Vaisala temperature dependence varies in magnitude and sign. Figure 3 shows the CO₂ mixing ratios and temperature dependence at the Montclair Elementary School site. Compared to the Picarro instrument, this site also demonstrates higher amplitude diurnal cycles (Fig. 3a), but these diurnal cycles are in phase with the reference instrument. Unlike EBMUD, the Montclair site exhibits a negative temperature dependence (Fig. 3c). Figure 3b shows the diurnal cycles at Montclair and the nearby node located at College Preparatory School (CPS). The comparison of these two sites suggests there may indeed be an amplification of the diurnal cycle at Montclair caused by a negative temperature dependence of the Vaisala instrument.

Figure 3 (a) CO₂ mixing ratios from May 2018 at Montclair and measured with a Picarro instrument at the Richmond Field Station supersite. (b) Montclair 2018 diurnal cycle compared with College Preparatory School. (c) Temperature dependence of CO₂ signal at Montclair.

2.4 Temperature correction method

The goal of our approach to accounting for the temperature dependence of the Vaisala instruments is to rely exclusively on the network itself and, if available, supplementary reference instruments, such as a Picarro, for derivation of correction factors to null sensor temperature dependence.

The method we developed builds on our method for accounting for drift in the instrument zero. To derive a temperature factor, we use hourly averaged [CO₂]_STP and node measurements of temperature (T). It is important to note that a major factor contributing to the temperature inside the node is whether the node is placed in the sun or shade. As a result, direct correlation with meteorological temperature measured outside the node is not strong. For a moving three-week window, at each hour (h), the lowest 10th percentile of [CO₂]_STP within ±1 ^∘C of T(h) is calculated. A running array of temperature-based 10th percentile data is created for both the Picarro supersite $\left({}^{\mathrm{Picarro}}[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{10}\mathit{\%}}^T\right)$ at the Richmond Field Station and each Vaisala instrument $\left({}^{\mathrm{Vaisala}}[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{10}\mathit{\%}}^T\right)$ using the temperature (T) of the Vaisala instrument. The Vaisala temperature is assumed to be the temperature that the instrument is responding to. $\left({}^{\mathrm{\Delta }}[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{10}\mathit{\%}}^T\right)$ is then calculated, where

$$\begin{array}{}\text{(3)}& {}^{\mathrm{\Delta }}\left[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{10}\mathit{\%}}^T{=}^{\mathrm{Vaisala}}\right[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{10}\mathit{\%}}^T{-}^{\mathrm{Picarro}}[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{10}\mathit{\%}}^T.\end{array}$$

A linear regression for ${}^{\mathrm{\Delta }}[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{10}\mathit{\%}}^T$ against T provides a slope (m_T) and intercept (b_T) for a moving three-week time window. We considered the possibility that the instrument response to temperature could be a zero shift and/or a change in the response to CO₂. We were able to achieve similar results assuming the temperature effect is entirely due to one or the other of these possibilities. As there is already substantial drift in the instrument zero, we proceed under the assumption that the effect can be entirely attributed to the temperature dependence of the instrument zero. The median of m_T (^medm_T)is calculated for the deployment period of the Vaisala sensor to determine the temperature-corrected offset and CO₂ mixing ratios of Vaisala CO₂ measurements, based on an additive error correction (Eqs. 4–5). When it is observed that either the offset bias, the temperature-dependent slope, or the time-dependent drift in the instrument zero shifts dramatically during a deployment period, the deployment is manually separated into different periods that are calibrated separately. The occurrence and magnitude of this varies between instruments (0–3 times during a two-year-long deployment), and is typically identified by routine checks for agreement between neighboring sensors. Shifts in the offset bias, the temperature-dependent slope, or the time-dependent drift appear as sudden or gradual offsets in mixing ratios measured by a sensor and its neighbors. Typical identified shifts in the offset bias, the temperature-dependent slope, or the time-dependent drift are on the order of 10 ppm, 0.5 ppm ^∘C⁻¹, and 0.1 ppm d⁻¹, respectively:

$$\begin{array}{}\text{(4)}& {\displaystyle }\left[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{offset}}^T{=}^{\mathrm{\Delta }}\right[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{10}\mathit{\%}}^T{-}^{\mathrm{med}}{m}_T\times T\text{(5)}& {\displaystyle }[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{corrected}}^T=[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{STP}}{-}^{\mathrm{med}}{m}_T\times T.\end{array}$$

An example calibration, demonstrating m_T and $[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{offset}}^T$ over time at EBMUD 2020, is shown in Fig. S1 in the Supplement. Following calculation of the temperature-corrected offset, the temporal drift slope and intercept of this corrected offset are calculated and corrected using the methods described above, resulting in the generation of the temperature- and drift-corrected CO₂ offset $\left([{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{offset}}^{T,\mathrm{drift}}\right)$.

The final temperature- and drift-corrected CO₂ $\left([{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{corrected}}^{T,\mathrm{drift}}\right)$ is then calculated as

$$\begin{array}{}\text{(6)}& [{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{corrected}}^{T,\mathrm{drift}}=[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{STP}}{-}^{\mathrm{med}}{m}_T\times T-[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{offset}}^{T,\mathrm{drift}}.\end{array}$$

Figure 4 (a) CO₂ mixing ratios measured by the Picarro instrument at the Richmond Field Station (blue solid), uncorrected CO₂ measured at Elsa Widenmann Elementary School (orange solid), and node temperature measured at Elsa (green dashed). (b) CO₂ mixing ratios at Elsa Widenmann Elementary School with no temperature correction (green), temperature correction applied (green) and measured with a Picarro instrument at the Richmond Field Station (blue). (c) Temperature dependence at Elsa Widenmann Elementary School of ${}^{\mathrm{\Delta }}[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{10}\mathit{\%}}^T$.

The majority of the BEACO₂N nodes examined demonstrated a strong linear relationship between $\left({}^{\mathrm{\Delta }}[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{10}\mathit{\%}}^T\right)$ and node temperature. However, the node at Elsa Widenmann Elementary School appeared to show a strong negative temperature dependence only on particularly warm days (Fig. 4a, c). The temperature dependence of $\left({}^{\mathrm{\Delta }}[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{10}\mathit{\%}}^T\right)$ for this node better fit a quadratic than a linear relationship. To account for nodes with a nonlinear temperature dependence, in cases where a quadratic fit improves the R² of the fit by more than 0.2, the $\left([{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{offset}}^T\right)$ and $\left([{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{corrected}}^{T,\mathrm{drift}}\right)$ are calculated via Eqs. (7)–(8):

$$\begin{array}{}\text{(7)}& {\displaystyle }\left[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{offset}}^T{=}^{\mathrm{\Delta }}\right[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{10}\mathit{\%}}^T{-}^{\mathrm{med}}{m}_T^{\mathrm{1}}\times T{-}^{\mathrm{med}}{m}_T^{\mathrm{2}}\times T^{\mathrm{2}}\text{(8)}& {\displaystyle }\begin{array}{rl}[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{corrected}}^{T,\mathrm{drift}}& =[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{STP}}{-}^{\mathrm{med}}{m}_T^{\mathrm{1}}\times T{-}^{\mathrm{med}}{m}_T^{\mathrm{2}}\times T^{\mathrm{2}}\\ & -[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{offset}}^{T,\mathrm{drift}},\end{array}\end{array}$$

where $m_T^{\mathrm{1}}$ and $m_T^{\mathrm{2}}$ are the first and second terms of the quadratic fit of ${}^{\mathrm{\Delta }}[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{10}\mathit{\%}}^T$ against T.

We attempted to determine a relationship between Vaisala sensor age and temperature-dependence slope, but m_T was only weakly correlated with sensor age (r≈0.3). We did, however find some evidence that older sensors had a larger likelihood of having a larger temperature dependence. For sensors less than 3 years since their initial deployment, 90 % had m_T<1 ppm ^∘C⁻¹ and 64 % had m_T<0.5 ppm ^∘C⁻¹. For sensors older than 3 years, 75 % had m_T<1 ppm ^∘C⁻¹. and 47 % had m_T<0.5 ppm ^∘C⁻¹.

3 Evaluation of calibration

Figures 5b, e, and 4c show the temperature dependence of ${}^{\mathrm{\Delta }}[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{10}\mathit{\%}}^T$ nodes located at EBMUD, Montclair, and Elsa Widenmann, respectively. Figures 5a, d, and 4b show a comparison of the data at EBMUD, Montclair, and Elsa, respectively, with and without adjustment for a temperature-dependent zero offset. With the application of the temperature correction, the magnitudes of the diurnal cycles are reduced and demonstrate much better agreement in amplitude and phase with the Picarro instrument. The resulting diurnal cycle at EMBUD shows a much more typical diurnal cycle for an urban site, with a maximum occurring at mid-morning (Fig. 5c). At Montclair, the magnitude of the diurnal cycle is reduced, reaching a maximum of ∼ 430 ppm CO₂ during the early morning, and a minimum of ∼ 412 ppm CO₂ during midday–a pattern much more aligned with the diurnal cycle exhibited at CPS (Figs. 3b, 5c).

Figure 5 CO₂ mixing ratios at (a) EBMUD and (d) Montclair with no temperature correction (orange), temperature correction applied (green), and measured with a Picarro instrument at the Richmond Field Station supersite (blue). Temperature dependence of ${}^{\mathrm{\Delta }}[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{10}\mathit{\%}}^T$ at (b) EBMUD and (e) Montclair. Diurnal cycle with and without temperature correction at (c) EBMUD and (f) Montclair.

Figure 6 Data from 2017 are in top panels, 2018 are in middle panels, 2020 are in bottom panels. (a) CO₂ mixing ratios at Laney College with no temperature correction (green), temperature correction applied (blue), and measured with a Picarro instrument at the Richmond Field Station supersite (2018 and 2020) or with the Exploratorium Buoy (2017). (b) Temperature dependence at Laney College of ${}^{\mathrm{\Delta }}[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{10}\mathit{\%}}^T$. (c) Laney College diurnal CO₂ cycle with (green) and without (orange) temperature correction.

Following confirmation of the effectiveness of the temperature correction method on the sensors deployed at EBMUD in 2020 (EBMUD 2020, hereafter sensors will be referred to with the name and site year) and Montclair 2018, we examined the temperature-corrected CO₂ data at the Laney College BEACO₂N site during the spring (March–June) of three different years when different Vaisala CarboCap GMP343 instruments were deployed. Given the hypothesis that the observed temperature dependence is due to temperature-dependent errors in the Vaisala CO₂ signal, a successful calibration should be sensor specific, rather than site specific. Figure 6b demonstrates the varying ${}^{\mathrm{\Delta }}[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{10}\mathit{\%}}^T$ temperature dependence during three different years with different instrument deployments. Each deployment has a distinct offset and slope of ${}^{\mathrm{\Delta }}[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{10}\mathit{\%}}^T$ vs. temperature. During all deployment years, the temperature correction results in better agreement between the reference instruments and the Vaisala data (e.g., 14 April 2020), while preserving local signals (e.g., 15 June 2018) (Fig. 6a). The correction is also effective for the data record before deployment of the Picarro reference instrument in August 2017, when the Exploratorium CO₂ Buoy, located in the San Francisco Bay, was used as a reference instrument (Fig. 6a). The correction of the CO₂ diurnal cycle at Laney College is most notable during 2017, although midday levels of CO₂ are reduced in the corrected data for 2018 and 2020 as well (Fig. 6c).

The temperature correction method was further validated by examining neighboring sites in two regions of the Bay Area during and before periods of high CO₂ during September 2020 northern California fires. The Richmond sites of Washington Elementary School, Nystrom Elementary School, Dejean Middle School, and Peres Elementary and the Vallejo sites of Beverly Hills Elementary School, Mare Island Health and Fitness Academy, Grace Patterson Elementary School, and Highland Elementary School were compared. The resulting temperature-dependent percent differences of CO₂ between adjacent sites are reduced to approximately 0 %–2 % from 1 %–5 % (Figs. S3, S6). Temperature corrections also result in better agreement in CO₂ mixing ratios between adjacent sites in Richmond (Figs. 7 and S2) and in Vallejo (Figs. S4, S5). The results were identical when a multiplicative correction term, rather than additive, was considered (e.g., if the temperature effect was assumed to be on the CO₂ signal magnitude rather than entirely on the instrument zero).

Figure 7 CO₂ mixing ratios during and before 2020 September wildfires at four adjacent sites in Richmond without (a) and with (b) temperature correction.

Comparison of nearest-neighbor sites

To assess the improvement in the network precision following application of the temperature-dependence correction, we combined observations from the entire Bay Area network using data from all of 2020. All sites with available data for more than one month of 2020 were included. Nearest neighbor pairs of each site were identified, where nearest neighbors to an individual site were considered as the closest BEACO₂N sites within a 2 km radius of the site. There are 53 unique nearest neighbor pairs.

For each nearest neighbor pair X and Y, an array of the fractional differences between sites was calculated as: $\left(\right[{\mathrm{CO}}_{\mathrm{2}}{]}_X-\left[{\mathrm{CO}}_{\mathrm{2}}{]}_Y\right)/[{\mathrm{CO}}_{\mathrm{2}}{]}_X$. This was done using both the measurements before and after correction for temperature-dependent instrument zero $\left(\left[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{corrected}}^{\mathrm{drift}}\text{ and }\right[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{corrected}}^{T,\mathrm{drift}}\right)$. Figure 8a and d show the fractional differences of each nearest neighbor pair as a histogram calculated using $[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{corrected}}^{\mathrm{drift}}$ and $[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{corrected}}^{T,\mathrm{drift}}$, respectively. Most nearest neighbor site pairs exhibit a distribution of fractional differences centered close to zero, with both positive and negative tails (Fig. 8a, d). The temperature correction results in a clear improvement of agreement between nearest neighboring sites, with the mean of the absolute value of the average fractional differences of all nearest neighbor pairs decreasing by a factor of 2 from 0.025 to 0.013. For $[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{corrected}}^{T,\mathrm{drift}}$, this represents an average difference of 6.5 ppm at [CO₂]=500 ppm. Figure 8b and e express the fractional differences of nearest neighbor pairs as a single distribution calculated using $[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{corrected}}^{\mathrm{drift}}$ and $[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{corrected}}^{T,\mathrm{drift}}$, respectively. Fit to a Lorentz distribution, the mean and scale parameter of the distribution of nearest neighbor pairs using $[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{corrected}}^{\mathrm{drift}}$ is 0.0026 and 0.014, respectively, without accounting for temperature dependence and there is a substantial narrowing of the distribution, resulting in a mean and scale parameter of 0.005 and 0.007, respectively, after accounting for the effect of a temperature-dependent offset.

Figure 8 Histogram of the fractional differences between nearest neighbors sites (a) without and (d) with the temperature correction applied. Different colors represent different pairs of neighboring sites. Histogram of the fractional differences between all aggregated nearest neighboring sites (b) without and (e) with the temperature correction applied fit to a Lorentz distribution. Network mean of the percent difference for each nearest neighbor pair averaged by 2 ^∘C bins (c) without and (f) with temperature correction applied.

Further analysis was performed to confirm that the temperature correction method eliminates any temperature-dependent disagreement between nearest neighboring sites. The nearest neighbor fractional differences of CO₂ data were separated into 2 ^∘C temperature bins. For each temperature bin, the absolute value of the mean fractional difference between each nearest neighbor pair, using either $[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{corrected}}^{\mathrm{drift}}$ or $[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{corrected}}^{T,\mathrm{drift}}$, was calculated. We then averaged the mean fractional difference in each temperature bin over all nearest neighbor pairs. A plot of the resulting network mean percent difference vs. temperature is shown in Fig. 8c and f, using $[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{corrected}}^{\mathrm{drift}}$ and $[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{corrected}}^{T,\mathrm{drift}}$ data, respectively. In the original data, the mean percent differences were greatest at both high and low temperatures. In the temperature-corrected data, there is no clear dependence of nearest neighbor mean percent differences on temperature. The mean percent difference at all temperatures is also reduced.

4 Assessment of the network error

Turner et al. (2016) suggested that a mismatch error of ∼ 1 ppm CO₂ would be compatible with relevant constraints on point, line, and area CO₂ sources of 147, 45, and 9 t C h⁻¹, respectively. Minimizing the network measurement error to close to 1 ppm is desirable, as at this measurement uncertainty, the error in emissions estimates from inverse modeling becomes dominated by model uncertainties (Turner et al., 2016). Assessing network error in the field is, however, a complex problem. We approach the problem by exploring differences between adjacent nodes, which should be an upper limit to the uncertainty. Although the site-to-site variation is strongly influenced by local emissions sources, there are also strong correlations with changes in urban-, synoptic-, and global-scale CO₂ signals that are spatially coherent across pairs of adjacent nodes. Variances between adjacent nodes are due to a combination of true site-specific signals and instrument biases. It is therefore difficult to know the minimum variance in adjacent nodes for a hypothetical “perfect” measurement. For nearest neighbor sites, the majority of the CO₂ signal should show near-zero difference, representing the background signal. In the observation record we would also expect moments when either site in a pair has a larger signal, driven by local emission sources and meteorology. Sites closer to the highway also typically have larger CO₂ signals (Shusterman et al., 2018). In the following section we describe a procedure for evaluating network error and summarize the improvements following inclusion of the temperature correction described above.

4.1 Site variance and correlation length scales

To evaluate the network error, a semivariogram of γ_nn vs. distance was constructed for $[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{corrected}}^{T,\mathrm{drift}}$ (Fig. 9). Using data from all sites with more than three days of available data during the summer of 2020, we calculated the semivariance between CO₂ measurements at each BEACO₂N node, i, and all other sites in the Bay Area network (Eq. 9)

$$\begin{array}{}\text{(9)}& {\mathit{\gamma }}_{nn}={\displaystyle \frac{{\sum }_j^N{\left([{\mathrm{CO}}_{\mathrm{2}}{]}_i-[{\mathrm{CO}}_{\mathrm{2}}{]}_j\right)}^{\mathrm{2}}}{\mathrm{2}N}}.\end{array}$$

Summer months were chosen because the average and diurnal variability of CO₂ mixing ratios are reduced, meaning that measured site variances are relatively more influenced by instrument error, rather than by “true” atmospheric variance, than in the winter. In Fig. 9 the square root of the semivariance is plotted against the distance separating the BEACO₂N nodes and fitted with an exponential model. The Picarro reference instrument at the Richmond Field Station was included in this analysis.

Figure 9 Semivariogram of BEACO₂N sites for data with temperature correction applied. Data are averaged by 0.1 km bins. Plot includes data from the Picarro instrument at Richmond Field Station.

Using the root semivariance as a correlation metric, in temperature-corrected data, the e-folding length scale for variation is 1.2 ± 0.3 km (1.7 ± 0.7 km using semivariance as a correlation metric, not shown), supporting the BEACO₂N hypothesis that 2 km node spacing in a dense network will capture important elements of local variability. The temperature-correction results in a maximum root semivariance of 5.5 ± 2 ppm (reduced from 8 ppm in the uncorrected data). Extrapolated to a distance of zero, the temperature correction method has a predicted root semivariance of 1.3 ± 0.9 ppm, representing the network error. This analysis suggests that the temperature correction method provides a meaningful reduction of network measurement uncertainties toward our desired ∼ 1 ppm network error.

Length scales for correlations (r²) between sites calculated by Shusterman et al. (2018) during the summer 2017 were larger than the 1.2 km length scale identified here for root semivariance (1.7 km for semivariance). To more directly compare, we also performed the method of Shusterman et al. (2018) on the temperature-corrected CO₂ data for the summer of 2020. We examined the correlation of CO₂ concentrations for every pairing of Bay Area sites during this period for all hours, during the day, and during the night (Fig. S7). The e-folding distance for the decay of r² correlation coefficients was 2.8 km for all times, 3.7 km during the day, and 2.8 km at night. This is in good agreement with the length scales of 2.9 km at all times, 3.6 km during the day, and 2.2 km at night found by Shusterman et al. (2018). The base-line correlation for sites separated by more than 20 km was found to be 0.46, larger than the correlation background of ∼0.3 of Shusterman et al. (2018). The temperature correction does not affect the characteristic length scale of BEACO₂N sites, but improves the overall base-line correlations and variances.

4.2 Contribution of instrument error to site variance

We can represent the network instrument error also by examining the sources contributing to the semivariance between nearest neighboring sites. The semivariance (γ_nn) of nearest neighboring sites can be expected to have contributions from both “true” variations in emissions and meteorology and erroneous differences caused by instrument error. “True” variations in emissions and meteorology are reflected in temporal changes in CO₂ concentrations due to emissions plumes and changes in wind speed and direction. Here we used temporal changes in CO₂ concentrations at a certain site as a proxy for “true” atmospheric variations in CO₂. To estimate the portion of the semivariance resulting from atmospheric phenomena, an analogous quantity for the hourly variations in CO₂ was calculated for each site according to Eq. (10):

$$\begin{array}{}\text{(10)}& {\mathit{\gamma }}_{hh}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}N}}\sum _{\mathrm{1}}^{N-\mathrm{1}}{\left([{\mathrm{CO}}_{\mathrm{2}}{]}_h-[{\mathrm{CO}}_{\mathrm{2}}{]}_{h+\mathrm{1}}\right)}^{\mathrm{2}},\end{array}$$

where N is the number of hours of data and [CO₂]_h , and [CO₂]_h+1 are the measured mixing ratios of CO₂ at each hour and 1 h later, respectively. The individual instrument error was then calculated as:

$$\begin{array}{}\text{(11)}& {\mathit{ϵ}}_{\mathrm{inst}}=\sqrt{{\mathit{\gamma }}_{nn}-{\mathit{\gamma }}_{hh}}.\end{array}$$

The resulting upper-bound instrument error from the median of individual instrument errors for the Bay Area network is 2.5 ± 0.5 ppm. (This estimate for nontemperature corrected data is 4.5 ± 0.9 ppm). We consider this an upper bound because hourly variations in the CO₂ signal reflects the atmospheric changes at an individual site, which may not match with the atmospheric changes at the nearest neighbor sites. Variations in emissions or wind velocity may result in larger “true” differences between a site and its nearest neighbor than are represented by the site's hourly variability.

To reduce the influence from short-term atmospheric variations, the network error was also estimated using an individual site's root mean squared error (RMSE_i) as a metric for “true” atmospheric variation (Eq. 12) and a “paired” RMSE (RMSE_paired) using the mean CO₂ signal of its nearest neighbor site ($\stackrel{\mathrm{‾}}{{}^{nn}\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$) as a measure of total variation (Eq. 13). The site error was then calculated according to Eq. 14:

$$\begin{array}{}\text{(12)}& {\displaystyle }{\mathrm{RMSE}}_i=\sqrt{{\displaystyle \frac{{\sum }_{h=\mathrm{1}}^N\left([{\mathrm{CO}}_{\mathrm{2}}{]}_h-\stackrel{\mathrm{‾}}{[{\mathrm{CO}}_{\mathrm{2}}{]}_i}\right)}{N}}}\text{(13)}& {\displaystyle }{\mathrm{RMSE}}_{\mathrm{paired}}=\sqrt{{\displaystyle \frac{{\sum }_{h=\mathrm{1}}^N\left([{\mathrm{CO}}_{\mathrm{2}}{]}_h-\stackrel{\mathrm{‾}}{{}^{nn}\left[{\mathrm{CO}}_{\mathrm{2}}\right]}\right)}{N}}}\text{(14)}& {\displaystyle }{\mathit{ϵ}}_{\mathrm{inst}}=\sqrt{{\mathrm{RMSE}}_{\mathrm{paired}}^{\mathrm{2}}-{\mathrm{RMSE}}_i^{\mathrm{2}}}.\end{array}$$

The resulting network instrument errors were between 0.5 and 4 ppm, with a median of 1.6 ± 0.4 ppm, in good agreement with the error calculated from the semivariogram fit. Based on these analyses, we estimate the network error of the Bay Area BEACO₂N network to be less than 1.6 ppm, close to our goal of 1 ppm network error.

5 Application to other city networks

The BEACO₂N network has recently been extended to several other cities, and will further expand to additional locations in coming years. Currently, locations where BEACO₂N nodes are deployed (in addition to the Bay Area) are Houston (19 nodes, network start November 2017), Glasgow in collaboration with the University of Strathclyde (> 20 nodes, network start May 2021), New York City (8 nodes, network start April 2018), and Los Angeles, in collaboration with the University of Southern California (12 nodes, network start May 2021). The goal of the network is to be self-calibrated, as not all locations at which the nodes will be deployed have a highly precise and frequently calibrated reference instrument. As such, an alternative method of obtaining a reference for the determination of drift, offset, and temperature dependence is needed.

We find that the network median [CO₂]_STP $\left([{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{STP}}^{\mathrm{med}}\right)$ can be used as a reference. To begin, we define the network median $\left([{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{STP}}^{\mathrm{med}}\right)$ as the median [CO₂]_STP of sites having a temperature-dependent slope (m_T) less than 1 ppm ^∘C⁻¹. $[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{STP}}^{\mathrm{med}}$ is used as a “reference site” from which temperature-based 10th percentile data $\left({}^{\mathrm{med}}[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{10}\mathit{\%}}^T\right)$ is calculated for the determination of ${}^{\mathrm{\Delta }}[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{10}\mathit{\%}}^T$:

$$\begin{array}{}\text{(15)}& {}^{\mathrm{\Delta }}\left[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{10}\mathit{\%}}^T{=}^{\mathrm{Vaisala}}\right[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{10}\mathit{\%}}^T{-}^{\mathrm{med}}[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{10}\mathit{\%}}^T.\end{array}$$

5.1 Bay Area tests

We observe good agreement between the Picarro reference instrument during 2020 and $[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{STP}}^{\mathrm{med}}$ (Fig. 10). The mean percent difference considering all 2020 data is 0.46 %, representing an accuracy error of 2 ppm at 420 ppm CO₂ (Fig. 10d). We also do not see evidence of a temperature-dependent offset between the Picarro reference instrument and $[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{STP}}^{\mathrm{med}}$.

Figure 10 (a) CO₂ mixing ratios measured by the Picarro instrument at the Richmond Field Station (blue) and the median CO₂ of all Bay Area nodes having a temperature-dependent slope less than an absolute value of 1 ppm ^∘C⁻¹ (orange). (b) The difference in the tenth percentile of CO₂ mixing ratios measured by the Picarro instrument and the network median plotted versus date and (c) versus temperature. (d) Histogram of the fractional differences between Picarro CO₂ mixing ratios and the network median. Data for (c) and (d) include all of 2020.

The precision of the Bay Area network is negligibly affected when the network median is used as the reference, with the mean of the absolute value of the average fractional differences of all nearest neighbor pairs equal to 0.015 ± 0.008 (compared to 0.013 ± 0.007 with the Picarro as reference) (Fig. S8). The resulting maximum root semivariance is 5.5 ± 2 ppm and extrapolated root semivariance at 0 km separation is 0.8 ± 0.9 ppm, respectively, approximately equal to the values calculated when the Picarro is used as a reference. The network accuracy is however, more appreciably altered. Figure 11 shows the fractional difference between $[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{corrected}}^{T,\mathrm{drift}}$ determined using the Picarro and $[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{STP}}^{\mathrm{med}}$ as a reference at each site. The resulting mean percent difference is 0.51 ± 0.02 %, representing a network accuracy error of 2 ppm at 420 ppm CO₂. This accuracy error is mainly driven by small differences in the offsets (2 ppm on average) and m_T (0.2 ppm ^∘C⁻¹ on average, see Supplement) between $[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{corrected}}^{T,\mathrm{drift}}$ calculated using the Picarro and $[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{STP}}^{\mathrm{med}}$ as a reference. These results suggest that the network precision can be expected to remain near 1 ppm CO₂ with the use of $[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{STP}}^{\mathrm{med}}$ as a reference, but additional accuracy error of 2 ppm may be introduced. The influence of a sea breeze in the Bay Area makes the median tenth percentile CO₂ measured by Bay Area nodes a regional background. Although the median tenth percentile of other inland sensor networks may not represent a regional background, it can be expected to represent the overall network regional average baseline.

Figure 11 (a) Histogram of the fractional differences between sites with temperature-corrected CO₂ calculated using the Picarro instrument and the Bay Area network median as a reference. Different colors represent different sites. The mean indicated is the average of the absolute values of each neighboring pair's mean fractional difference. (b) Histogram of the aggregated fractional differences between sites with temperature-corrected CO₂ calculated using the Picarro instrument and the Bay Area network median as a reference. The mean and error indicated are the mean and 95 % confidence interval of the distribution.

Analysis of the Bay Area network was performed on the 36 nodes with sufficient data availability for 2020. However, the newly established networks have fewer nodes than in the Bay Area. To use $[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{STP}}^{\mathrm{med}}$ as a reference, we must have sufficient nodes from which to calculate the network median. To evaluate this, for n=1–26, a random subset of n Bay Area nodes was selected 100 times. For each of the 100 random subsets of n nodes, the mean fraction difference was calculated between the network median CO₂ and the median calculated using the subset. The average and standard error of the 100 mean fraction differences was then calculated. The results of this analysis are presented in Fig. 12. We suggest that a minimum of 7 nodes with m_T less than 1 ppm ^∘C⁻¹ is required for the accuracy error to be lower than 2 %. For less than 1 % error, at least 12 nodes are required.

Figure 12 Fractional difference between the Bay Area network median calculated from all Bay Area sites and the network median calculated from a subset of between one and 26 nodes. A random subset of n=1–26 nodes were selected to calculate the mean fractional difference between the network median CO₂ and the median calculated using the subset. This was repeated 100 times for each of n=1–26 nodes. The reported fractional difference and error bars are the average and 95 % confidence interval of the mean fractional difference from the 100 random samples.

5.2 Houston

Data from the Houston network were subsequently calibrated using $[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{STP}}^{\mathrm{med}}$ as a reference for determining temperature dependence, drift, and offset. Temperature dependence calibration of each site in the Houston network was performed twice. All sites were first included in $[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{STP}}^{\mathrm{med}}$ and sites with m_T greater than 1 ppm ^∘C⁻¹ were identified. These sites were then excluded from $[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{STP}}^{\mathrm{med}}$ and each site was recalibrated. Histograms of the fraction differences between nearest neighbor sites are shown in Fig. 13. The average mean percent difference between nearest neighbors was 2 ± 1 %. Though considerably larger than the differences between nearest neighbors in the Bay Area network, it is not immediately clear whether this difference is caused by greater precision error in Houston or differing meteorology and CO₂ sources that cause greater differences between CO₂ mixing ratios at adjacent sites. We attempted to perform a similar instrumental error analysis, but there are currently insufficient overlapping CO₂ data in Houston for uncertainty analysis. However, we do not have reason to expect the instrument errors would be any larger in the Houston network.

Figure 13 (a) Histogram of the fractional differences between nearest neighbors sites in the Houston network with the temperature correction applied using the network median as a reference. Different colors represent different pairs of neighboring sites. The mean indicated is the average of the absolute values of each neighboring pair’s mean fractional difference. The error is the associated 95 % confidence interval. (b) Histogram of the fractional differences between all aggregated nearest neighbors sites with the temperature correction applied. The mean and error indicated are the mean and 95 % confidence interval, respectively, of the distribution.

6 Conclusions

We have assessed the accuracy of the BEACO₂N network following in situ calibration of the temperature-dependence in Vaisala CO₂ sensors. We report meaningful reductions in network uncertainties following application of a temperature-dependence correction, and a resulting network instrument error of 1.6 ppm CO₂ or less.

A method for correcting Vaisala instrument temperature dependence in BEACO₂N has been established and evaluated using sites across the San Francisco Bay Area network. The method corrects observations from individual instruments so that they exhibit a temperature dependence in their lowest temperature-based 10th percentile of CO₂ data that is equivalent to that of a reference site, thus correcting erroneous instrument temperature dependence while preserving true diurnal cycles and local signals. This field calibration of temperature dependence can be entirely internal to the network and does not necessarily require a reference instrument, although the addition of a reference instrument provides greater network accuracy. The implementation of the temperature correction method produces more reasonable diurnal cycles, diurnal cycles that are maintained for sites influenced by similar emissions sources, and better agreement between adjacent sites. We additionally describe methods for characterizing network scale uncertainties and site-to-site biases. The average variation between adjacent sites was found to be 1.3 % following implementation of temperature correction (compared to 2.5 % prior to the correction). The temperature correction greatly improves the precision of CO₂ measurements in the BEACO₂N network.

We show that the network precision can be maintained at 1.3 % even in locations without a high-cost reference instrument, using the network median as a reference, provided that there are at least 12 sites with small temperature dependencies. This has important implications for the expansion of BEACO₂N to additional cities globally, as well as for other dense low-cost CO₂ networks. However, without a reference instrument, the network accuracy error increases relative to a network that utilizes a reference instrument by $\sim \pm$2 ppm.