Revised methodology for CO₂ and CH₄ measurements at remote sites using a working standard-gas-saving system

We have revised a calculation method of mole fractions and uncertainties for in situ CO₂ and CH₄ measurements with a working standard-gas-saving system. It uses on-site compressed air to track the baseline drift of sensors. The Japan–Russia Siberian Tall Tower Inland Observation Network (JR-STATION) is made up of this system, which was installed across nine different sites in Siberia. The system acquires semi-continuous data by alternating between sampling air from multiple altitudes through switched flow paths and recording several minutes of averaged data for each altitude. We estimated the sensor repeatability (u_r) based on the measurement of on-site compressed air. The u_r for CO₂ and CH₄ was mostly around 0.05 ppm and below 5 ppb, respectively. The combined standard uncertainties (u_c(x)) of time-averaged ambient air measurements were sometimes higher than the u_r for each period because the data included atmospheric variability during the measurement period of several minutes. Data users should consider the difference between the u_r and u_c(x) to select optimal data, depending on their focusing spatial scale. The CO₂ and CH₄ data measured with a non-dispersive infrared (NDIR) analyzer and a tin dioxide sensor (TOS) exhibited good agreement with those measured by cavity ring-down spectroscopy (CRDS).

1 Introduction

It is known that accurate measurements of greenhouse gas mole fractions require the analyzers to be calibrated against a set of standard gas mixtures. At least one of them (target) should be used hourly to track a non-dispersive infrared (NDIR) analyzer's baseline drift (Andrews et al., 2014). Delivering high-pressure cylinders to remote sites is a significant issue for long-term atmospheric monitoring. Thus, to reduce the consumption of gas, Watai et al. (2010) developed a system that utilizes on-site air as a sub-working standard gas (SWS-gas) to track the baseline drift of the NDIR sensors. The National Institute for Environmental Studies (NIES) installed this system at a remote tower site at Berezorechka (56°08^′45^′′ N, 84°19^′55^′′ E) in western Siberia in 2001 to measure CO₂ mole fractions semi-continuously. After this, in central Siberia, Winderlich et al. (2010) developed a measuring system without dehumidification using a cavity ring-down spectroscopy (CRDS) analyzer to reduce the frequency of cylinder replacement. The CRDS instrument is a more stable device, and a calibration frequency of every 2 weeks to every month is recommended (ICOS RI, 2020).

Concerning CH₄ measurement, a commonly used gas chromatograph with a flame ionization detector (GC/FID) requires hydrogen and carrier gases. It also needs significant power consumption. However, electric power is often restricted at remote sites. Thus, Suto and Inoue (2010) modified a tin dioxide sensor (TOS), which is widely used to detect natural gas leaks, to be able to measure CH₄ in the atmosphere. The developed TOS unit does not need hydrogen and carrier gases. The nominal power consumption for the unit, consisting of the TOS, a temperature stabilizer mechanism, and electronic circuits for the sensor and heater, is less than 10 W.

We added the TOS unit to the system at the tower site in western Siberia, then expanded the tower observation network (Sasakawa et al., 2010, 2012, 2013). The tower network named JR-STATION now consists of six tower sites in western Siberia. Recently, we added CRDS analyzers at Karasevoe (58°14^′44^′′ N, 82°25^′28^′′ E) in 2015 (Picarro G2401) and at Demyanskoe (59°47^′29^′′ N, 70°52^′16^′′ E) and Noyabrsk (63°25^′45^′′ N, 75°46^′48^′′ E) in 2016 (Picarro G2301) to improve the robustness of the measurements.

The system follows the operational method conceived under the instrumentation constraints a quarter of a century ago (around 2000) and in remote areas with limited infrastructure (Watai et al., 2010). However, Watai et al. (2010) did not present an optimal calculation method for the measurement sequence of this system (nor did the TOS calculation method). Nor did they calculate uncertainties, especially as has been recommended in the GAW (Global Atmosphere Watch) report (2020) in recent years.

Thus, we have updated the method for calculating CO₂ and CH₄ mole fractions to derive their uncertainty for each dataset simultaneously. Here, we describe the details of the modified measurement system and the calculation method. Furthermore, some sites have installed CRDS instruments, which have been widely used for greenhouse gas observations (Yver Kwok et al., 2015) and allow partial comparisons with conventional sensors, so the recalculated data were compared with the CRDS data to see how well they agree.

2 Method

2.1 Measurement system

Ambient air was taken from air sample inlets at two different heights (four at Berezorechka) on television and radio-relay communication towers (Table 1). Each sample inlet was mounted several meters away from the tower at the end of an extension arm. The air from the inlets was pulled into the measurement system through the sampling lines (6 mm outer diameter Decabon tube). The measurement system was housed in a freight container insulated to reduce temperature variation. Two thermometers were mounted inside the container, one near the ceiling and the other near the floor. According to the upper thermometer, the room temperature in the container was kept above 15 °C and the temperature difference in the 12 h calibration interval was kept below 3 °C on average during the year. Since the introduction of CRDS, a simple cooler was installed to prevent the temperature inside the container from rising too high during the summer months due to the heat generated by CRDS. A schematic diagram of the measurement system is shown in Fig. 1. The measurement system consists of a pump unit, a selector unit, and an analyzer unit. The pump unit was located upstream of the selector and analyzer unit to keep the downstream pressure higher than ambient, which reduced the likelihood of bias in measurements due to any leak from many connections in the system. Two diaphragm pumps (model N86KNE, KNF, Germany) delivered the sample air into the system. The sampling lines were flushed continuously with a flow rate of about 7 standard liters per minute, and excess air was vented through the back-pressure valve (BPV in the pump unit). Then the air was dried by an adiabatic expansion in a glass water trap (WT in the pump unit) that was purged every hour via an NC solenoid valve, which was opened twice for 10 s to remove the condensed water. The sample air was also dried using a semipermeable membrane dryer (PD-625-24SS, Perma Pure, USA) (Nafion in the selector unit). The semipermeable membrane dryer removed water vapor from the pressurized inner tube to an outer tube where the split gas flowed (split sample method). The air from the upper- and lower-level inlet, the three working standard gases (WS-gases), and the sub-working standard gas (SWS-gas) were selected through a six-port valve manifold. While the WS-gases or the SWS-gas flowed into the analyzer unit, the sample air was exhausted at the six-port valve. In the analyzer unit, the sampled air was extra-dried with magnesium perchlorate, which was fed into a stainless-steel tube with a dimension of 2 cm in inner diameter and 10 cm in length (Mg(ClO₄)₂ in the analyzer unit). There were two tubes, and the flow path of the air switched from one to the other every month. The used magnesium perchlorate was replaced before the next run. After being dried with the magnesium perchlorate, the air retained its dew point at around −50 °C (39 ppm). The dehumidified air was then introduced into an NDIR analyzer (LI-820, LI-COR, USA; LI-7000 was used until September 2008 at BRZ) at a constant flow rate of 35 standard cubic centimeters per minute (sccm) set by a mass flow controller (SEC-E40, STEC, Japan). The CO₂ mole fraction was defined as the mole fraction in the dry air, and water vapor correction was not adopted. After passing through the NDIR analyzer, the air flowed into the TOS unit. A chemical desiccant made of P₂O₅ was installed in front of the TOS because it is necessary to keep water vapor below 10 ppm in the sample air for this type of sensor. The sensor was designed to operate in areas lacking the sufficient infrastructure to sustain a conventional measurement system, such as a significant power source, carrier gas supply, and temperature-stabilized environment. The sensor has been verified against a gas chromatograph equipped with a flame ionization detector (Suto and Inoue, 2010). We additionally installed the CRDS (Picarro Inc.) analyzer at Karasevoe (G2401) in 2015 and at Demyanskoe and Noyabrsk (G2301) in 2016 to improve the system. The sampled air was split after leaving the six-port valve, then fed into the CRDS instrument at a constant flow rate of 35 sccm set by a mass flow controller (SEC-E40, STEC, Japan) through a semipermeable membrane dryer (model MD-050-72S-1, Perma Pure, USA). To protect the cavity of the CRDS instrument from an inflow of the dissolved chemical desiccant (Mg(ClO₄)₂ or P₂O₅) in the accidental case of a broken pump, we equipped the CRDS instrument with two poppet check valves (PCV in the analyzer unit). When the pumps in the pump unit stop and only the CRDS pump is running, the flow stops at the PCV upstream of P₂O₅, and the increased suction pressure allows air in the container to enter from the PCV in front of the Nafion dryer. The data from this process have been deleted.

Figure 1 Schematic diagram of the tower observation system.

Three WS-gases (STD1, STD2, STD3) were prepared from pure CO₂ and CH₄ (G1 grade, Japan Fine Products Corp. (JFP), Japan) diluted with purified air (G1 grade, JFP), and their mole fractions were determined against the NIES 09 CO₂ scale (Machida et al., 2011) and NIES 94 CH₄ scale. Each scale was established by a series of standard gases prepared by the gravimetric method. Since the pure CO₂ gas is derived from burned petroleum, the isotopic CO₂ composition of the gases is lighter than atmospheric CO₂. When the NDIR analyzer is calibrated against CO₂ standards with a lighter-than-atmospheric-CO₂ isotopic composition, the NDIR analyzer measures a lower CO₂ mole fraction in sample air with a known CO₂ mole fraction. The error in the apparent NDIR CO₂ mole fraction depends on its individual sensitivity to the optical filter property. Tohjima et al. (2009) reported that the errors for the three NDIR analyzers range from −0.04 to −0.08 ppm. Compared to the WMO-CO2-X2007 scale, the NIES 09 CO₂ scale is consistent within 0.1 ppm (Round Robin 5 and 6 Comparison Experiment). Since there have yet to be published results with WMO-CO2-X2019, scale conversion could be the linear function shown in Hall et al. (2021) (X2019 = 1.00079 × NIES 09 − 0.142 ppm). The NIES 94 CH₄ scale ranges from 3.0 to 5.5 ppb higher than the WMO-CH4-X2004A scale (Round Robin 5 and 6 Comparison Experiment).

Table 1 Main features of tall towers in the network used for continuous long-term atmospheric CO₂ and CH₄ measurements over Siberia.

^* Approximate estimates from Google Earth.

A frequent calibration with WS-gases within 1–2 h is necessary to conduct precise measurements of CO₂ and CH₄ because the output of the NDIR analyzer or the TOS could vary depending on the environment (e.g., atmospheric pressure) in 1–2 h. But if the calibration were done at this frequency, standard gases would be consumed in less than a year. Because delivering WS-gases to remote sites is a significant issue, we utilized on-site compressed air as an SWS-gas to track the sensors' baseline drift, which reduced the consumption of the three standard gases. The on-site compressed air (“On-site gas A/B” in Fig. 1) was analyzed every hour, and the WS-gases were measured every 12 h to calibrate the sensors (details of the sequence are shown below). The measurement protocol adhered to the procedure established by Watai et al. (2010) for NDIR analysis using this system. However, at sites where CRDS was installed, the WS-gas measurement interval was extended to 48 h to prolong the longevity of the WS-gas. An aluminum cylinder (0.048 m³) for SWS-gas was automatically exchanged when the inner pressure decreased below 0.1 MPa; then soon air from the highest inlet was compressed by a pump (LOA-P103-NO, GAST, USA) into the cylinder for about 5 h to approximately 0.35 MPa, after having been passed through a similar triple dehumidification path as the sampled air (a stainless-steel water trap, a semipermeable membrane dryer (SWF-M06–400, AGC, Japan), and magnesium perchlorate). It was preserved for approximately 1 week (3 d with CRDS) for usage until the inner pressure in one used for measurements decreased below 0.1 MPa (Table 2). Schibig et al. (2018) reported that the CO₂ mole fraction in a 29.5 L aluminum cylinder increases by 0.090 ± 0.009 µmol mol⁻¹ when dropping from 150 to 1 bar but also note that this change is smaller if larger cylinders are used. Given the cylinder size and filling pressure in this system, mole fraction changes within the cylinder are considered negligible. The variations in SWS mole fraction with CRDS between WS-gases (48 h) were in fact very stable regardless of the SWS mole fraction range (Figs. S1–S2 in the Supplement).

Table 2 SWS-gas A/B measurement and filling sequence.

Air temperature and relative humidity were measured at both heights on the tower using commercial sensors (HMP45D, Vaisala, Finland). A wind monitor (model 81000, R. M. Young, USA) determined wind direction and speed at the higher inlet. Solar radiation was measured by a pyranometer (CM3, Kipp & Zonen, the Netherlands) and precipitation by a tipping-bucket rain gauge (model 52202, R. M. Young, USA) on the top of the container laboratory. The analysis operation and data logging were performed by a measurement and control system (CR10X datalogger, CAMPBELL, USA). Stored data were retrieved once a month when both a system check and replacement of consumables (e.g., chemical desiccants) took place.

2.2 Measurement sequence

To be able to measure air at two heights, the air-sampling flow path was rotated every 20 min with the six-port valve in the selector unit; that is, the higher inlet was sampled from hour 00 to hour 20, the lower inlet from hour 20 to hour 40, and the SWS-gas from hour 40 to (hour +1) 00. During the first 17 min of each 20 min sampling interval, the system is flushed to equilibrate the air sample after switching. The final 3 min readouts were averaged and reported as the representative output data for the applicable 1 h period. Measurement frequency was 3 s; thus, only the average and standard deviation (SD) of 60 readouts in voltage were stored in the CR10X. This was to minimize the data size for the limited storage capacity. The timestamp was the end time of every 20 min measurement interval. The raw data collected with the CRDS analyzer were stored in the CRDS instrument's hard disk and processed after downloading in our laboratory.

Figure 2 shows the schematic measurement sequence for half a day. In Fig. 2, we defined when the SWS-gas was measured just before an arbitrary series of WS-gas measurements as t₀. Then we numbered the time of the following measurements in turn. We also defined the series of standard gas measurements at the beginning of the 12 h as “B” and at the end as “E”.

Figure 2 Measurement sequence for a half-day between subsequent measurements of WS-gases.

2.3 Quality check of the standard gas measurements

First, we checked the relationship among three standard gas measurements. We calculated the differences (Δ^B(t_i), Δ^E(t_j)) between the measured output voltages of the standard gases (V_std(t_i), V_std(t_j)) and the estimated one of the SWS-gas at the time of the standard gas measurement (Fig. 2). Here i = 1, 2, 3 and j = 37, 38, 39. The output value of the SWS-gas was interpolated by time using the closest output of the SWS-gas before and after the series of standard gas measurements. Thus, these values and their variances are expressed as follows.

$$\begin{array}{}\text{(1)}& {\mathrm{\Delta }}^{\mathrm{B}}\left(t_i\right)=V_{\mathrm{std}}\left(t_i\right)-\left({\displaystyle \frac{\mathrm{6}-i}{\mathrm{6}}}\cdot V_{\mathrm{sub}}\left(t_{\mathrm{0}}\right)+{\displaystyle \frac{i}{\mathrm{6}}}\cdot V_{\mathrm{sub}}\left(t_{\mathrm{6}}\right)\right)\end{array}$$

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\left({\mathit{\sigma }}^{\mathrm{B}}\left(t_i\right)\right)}^{\mathrm{2}}& ={\left({\mathit{\sigma }}_{\mathrm{std}}\left(t_i\right)\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\mathrm{6}-i}{\mathrm{6}}}\cdot {\mathit{\sigma }}_{\mathrm{sub}}\left(t_{\mathrm{0}}\right)\right)}^{\mathrm{2}}\\ & +{\left({\displaystyle \frac{i}{\mathrm{6}}}\cdot {\mathit{\sigma }}_{\mathrm{sub}}\left(t_{\mathrm{6}}\right)\right)}^{\mathrm{2}}\end{array}\end{array}$$

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\mathrm{\Delta }}^{\mathrm{E}}\left(t_j\right)& =V_{\mathrm{std}}\left(t_j\right)\\ & -\left({\displaystyle \frac{\mathrm{42}-j}{\mathrm{6}}}\cdot V_{\mathrm{sub}}\left(t_{\mathrm{36}}\right)+{\displaystyle \frac{j-\mathrm{36}}{\mathrm{6}}}\cdot V_{\mathrm{sub}}\left(t_{\mathrm{42}}\right)\right)\end{array}\end{array}$$

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\left({\mathit{\sigma }}^{\mathrm{E}}\left(t_j\right)\right)}^{\mathrm{2}}& ={\left({\mathit{\sigma }}_{\mathrm{std}}\left(t_j\right)\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\mathrm{42}-j}{\mathrm{6}}}\cdot {\mathit{\sigma }}_{\mathrm{sub}}\left(t_{\mathrm{36}}\right)\right)}^{\mathrm{2}}\\ & +{\left({\displaystyle \frac{j-\mathrm{36}}{\mathrm{6}}}\cdot {\mathit{\sigma }}_{\mathrm{sub}}\left(t_{\mathrm{42}}\right)\right)}^{\mathrm{2}}\end{array}\end{array}$$

We estimated the output of STD1 at t₂ by adding Δ^B(t₁) to the estimated one of the SWS-gas at t₂. We also evaluated the output of STD3 at t₂ by adding Δ^B(t₃) to the estimated one of the SWS-gas at t₂. The same estimation was done at t₃₈. We then made a linear calibration line with the output of STD2 and the estimated outputs of STD1 and STD3. Only sets of the three standard gas measurements whose coefficients of determination were higher than 0.999 for CO₂ and 0.99 for CH₄ were adopted for the following calculation.

The difference in output (voltage) between Δ^B and Δ^E for each standard gas was defined as follows.

$$\begin{array}{}\text{(5)}& {\mathit{\delta }}_{\left(\mathrm{1},\mathrm{37}\right)}={\mathrm{\Delta }}^{\mathrm{E}}\left(t_{\mathrm{37}}\right)-{\mathrm{\Delta }}^{\mathrm{B}}\left(t_{\mathrm{1}}\right)\end{array}$$

$$\begin{array}{}\text{(6)}& {\mathit{\delta }}_{\left(\mathrm{2},\mathrm{38}\right)}={\mathrm{\Delta }}^{\mathrm{E}}\left(t_{\mathrm{38}}\right)-{\mathrm{\Delta }}^{\mathrm{B}}\left(t_{\mathrm{2}}\right)\end{array}$$

$$\begin{array}{}\text{(7)}& {\mathit{\delta }}_{\left(\mathrm{3},\mathrm{39}\right)}={\mathrm{\Delta }}^{\mathrm{E}}\left(t_{\mathrm{39}}\right)-{\mathrm{\Delta }}^{\mathrm{B}}\left(t_{\mathrm{3}}\right)\end{array}$$

The δ must be small unless the system is unstable, e.g., when the sensitivity of the sensors changes considerably for some reason. To exclude the data obtained during system malfunction, we determined a threshold for δ by converting it into mole fraction (< 5.0 ppm for CO₂, < 50 ppb for CH₄). Data showing values over the threshold were excluded from the calculation. The difference in CO₂ mole fraction between sides B and E was calculated as follows:

$$\begin{array}{}\text{(8)}& \stackrel{\mathrm{‾}}{{\mathit{\delta }}_{(i,j)}^{\mathrm{B}}}=\left|{\mathit{\delta }}_{(i,j)}/S^{\mathrm{B}}\right|,\end{array}$$

$$\begin{array}{}\text{(9)}& \stackrel{\mathrm{‾}}{{\mathit{\delta }}_{(i,j)}^{\mathrm{E}}}=\left|{\mathit{\delta }}_{(i,j)}/S^{\mathrm{E}}\right|,\end{array}$$

where S^B and S^E are the slopes of the linear regression line at sides B and E. Because the x axis of the calibration line for CH₄ is the logarithm of the mole fraction, the difference in CH₄ mole fractions was calculated as

$$\begin{array}{}\text{(10)}& \stackrel{\mathrm{‾}}{{\mathit{\delta }}_{(i,j)}^{\mathrm{B}}}=C_i\cdot \left|e^{\frac{{\mathit{\delta }}_{(i,j)}}{S^{\mathrm{B}}}}-\mathrm{1}\right|,\end{array}$$

$$\begin{array}{}\text{(11)}& \stackrel{\mathrm{‾}}{{\mathit{\delta }}_{(i,j)}^{\mathrm{E}}}=C_i\cdot \left|e^{-\frac{{\mathit{\delta }}_{(i,j)}}{S^{\mathrm{E}}}}-\mathrm{1}\right|,\end{array}$$

where C_i is the mole fraction of the standard gas.

2.4 Calculation of the sample mole fraction and the combined standard uncertainty

The analysis precision for this system under laboratory conditions was uniformly estimated as 0.3 ppm for CO₂ (Watai et al., 2010). Concerning CH₄ precision, Sasakawa et al. (2010) estimated it as 3.0 ppb based on the result of Suto and Inoue (2010). However, the experiment condition from Suto and Inoue (2010) was different from the gas-saving system. Instead, they connected only the WS-gases to the TOS, then reported the SD of repeated measurements. The CH₄ analysis precision for this system could thus be more significant than 3.0 ppb. Furthermore, the sensitivity and stability of the sensor could differ depending on the individual sensor and the condition of the individual system. We have thus updated the method for calculating the CO₂ and CH₄ mole fractions to derive their combined standard uncertainty for all data simultaneously.

2.4.1 Estimation of the output values of working standard gases and their SD at the time of the air sample measurements

We estimated the outputs in voltage of three standard gases at each measurement time of the sample air by interpolating the outputs of the three WS-gases depending on the difference of the outputs in voltage of the SWS-gas only when both standard gas measurements satisfied the criteria described in Sect. 2.3. Depending on the time difference between the targeted sample (t_k; k = 4, 5, 7, 8, … ,34, 35) and standard gases at both sides B and E ((t_i,t_j); (i,j) = $(\mathrm{1},\mathrm{37}),(\mathrm{2},\mathrm{38}),(\mathrm{3},\mathrm{39})$), the representative value $\left({\widehat{V}}_{\mathrm{std}\left(t_i,t_j\right)}^{\mathrm{BE}}\left(t_k\right)\right)$ and its variance $\left({\left({\widehat{\mathit{\sigma }}}_{\mathrm{std}\left(t_i,t_j\right)}^{\mathrm{BE}}\left(t_k\right)\right)}^{\mathrm{2}}\right)$ were estimated as follows (Fig. 3).

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& {\widehat{V}}_{\mathrm{std}\left(t_i,t_j\right)}^{\mathrm{BE}}\left(t_k\right)={\displaystyle \frac{a}{a+b}}\cdot {\widehat{V}}_{\mathrm{std}\left(t_i\right)}^{\mathrm{B}}\left(t_k\right)+{\displaystyle \frac{b}{a+b}}\cdot {\widehat{V}}_{\mathrm{std}\left(t_j\right)}^{\mathrm{E}}\left(t_k\right)\\ & ={\widehat{V}}_{\mathrm{sub}}\left(t_k\right)+{\displaystyle \frac{\mathrm{1}}{a+b}}\cdot \left\{a\cdot {\mathrm{\Delta }}^{\mathrm{B}}\left(t_i\right)+b\cdot {\mathrm{\Delta }}^{\mathrm{E}}\left(t_j\right)\right\}\end{array}\end{array}$$

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& {\left({\widehat{\mathit{\sigma }}}_{\mathrm{std}\left(t_i,t_j\right)}^{\mathrm{BE}}\left(t_k\right)\right)}^{\mathrm{2}}={\left({\widehat{\mathit{\sigma }}}_{\mathrm{sub}}\left(t_k\right)\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{a}{a+b}}\cdot {\mathit{\sigma }}_{\mathrm{std}}\left(t_i\right)\right)}^{\mathrm{2}}\\ & +{\left({\displaystyle \frac{a}{a+b}}\cdot {\displaystyle \frac{\mathrm{6}-i}{\mathrm{6}}}\cdot {\mathit{\sigma }}_{\mathrm{sub}}\left(t_{\mathrm{0}}\right)\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{a}{a+b}}\cdot {\displaystyle \frac{i}{\mathrm{6}}}\cdot {\mathit{\sigma }}_{\mathrm{sub}}\left(t_{\mathrm{6}}\right)\right)}^{\mathrm{2}}\\ & +{\left({\displaystyle \frac{b}{a+b}}\cdot {\mathit{\sigma }}_{\mathrm{std}}\left(t_j\right)\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{b}{a+b}}\cdot {\displaystyle \frac{\mathrm{42}-j}{\mathrm{6}}}\cdot {\mathit{\sigma }}_{\mathrm{sub}}\left(t_{\mathrm{36}}\right)\right)}^{\mathrm{2}}\\ & +{\left({\displaystyle \frac{b}{a+b}}\cdot {\displaystyle \frac{j-\mathrm{36}}{\mathrm{6}}}\cdot {\mathit{\sigma }}_{\mathrm{sub}}\left(t_{\mathrm{42}}\right)\right)}^{\mathrm{2}}\end{array}\end{array}$$

Here, a : b = (t_j−t_k) : (t_k−t_i). The hat ( $\widehat{}$) means the estimated value. ${\widehat{V}}_{\mathrm{sub}}\left(t_k\right)$ was calculated by interpolating the output of the SWS-gas value nearest to the targeted sample as follows.

$$\left\{\text{spl}\mathrm{1}\mathrm{|}k=\mathrm{4}\right\}$$

$$\begin{array}{}\text{(14)}& {\widehat{V}}_{\mathrm{sub}}\left(t_{\mathrm{4}}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\cdot V_{\mathrm{sub}}\left(t_{\mathrm{0}}\right)+{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}\cdot V_{\mathrm{sub}}\left(t_{\mathrm{6}}\right)\end{array}$$

$$\left\{\text{spl}\mathrm{1}\mathrm{|}k=\mathrm{7},\mathrm{10},\mathrm{13},\mathrm{\dots },\mathrm{34}\right\}$$

$$\begin{array}{}\text{(15)}& {\widehat{V}}_{\mathrm{sub}}\left(t_k\right)={\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}\cdot V_{\mathrm{sub}}\left(t_{k-\mathrm{1}}\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\cdot V_{\mathrm{sub}}\left(t_{k+\mathrm{2}}\right)\end{array}$$

$$\left\{\text{spl}\mathrm{2}\mathrm{|}k=\mathrm{5}\right\}$$

$$\begin{array}{}\text{(16)}& {\widehat{V}}_{\mathrm{sub}}\left(t_{\mathrm{5}}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{6}}}\cdot V_{\mathrm{sub}}\left(t_{\mathrm{0}}\right)+{\displaystyle \frac{\mathrm{5}}{\mathrm{6}}}\cdot V_{\mathrm{sub}}\left(t_{\mathrm{6}}\right)\end{array}$$

$$\left\{\text{spl}\mathrm{2}\mathrm{|}k=\mathrm{8},\mathrm{11},\mathrm{14},\mathrm{\dots },\mathrm{35}\right\}$$

$$\begin{array}{}\text{(17)}& {\widehat{V}}_{\mathrm{sub}}\left(t_k\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\cdot V_{\mathrm{sub}}\left(t_{k-\mathrm{2}}\right)+{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}\cdot V_{\mathrm{sub}}\left(t_{k+\mathrm{1}}\right)\end{array}$$

Here $V_{\mathrm{sub}}\left(t_l\right)\left\{\text{sub}\mathrm{|}l=\mathrm{0},\mathrm{6},\mathrm{9},\mathrm{12},\mathrm{\dots },\mathrm{36}\right\}$ is the measured value. In the following, a calculation example of ${\widehat{V}}_{\mathrm{std}\left(t_i,t_j\right)}^{\mathrm{BE}}\left(t_k\right)$ and the variance for the case $\left\{\text{spl}\mathrm{1}\mathrm{|}k=\mathrm{4},\left(i,j\right)\right\}$ are given without any estimated value.

$$\left\{\text{spl}\mathrm{1}\mathrm{|}k=\mathrm{4},\left(i,j\right)=\left(\mathrm{1},\mathrm{37}\right),\left(\mathrm{2},\mathrm{38}\right),\left(\mathrm{3},\mathrm{39}\right)\right\}$$

$$\begin{array}{}\text{(18)}& \begin{array}{rl}& {\widehat{V}}_{\mathrm{std}\left(t_i,t_j\right)}^{\mathrm{BE}}\left(t_{\mathrm{4}}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\cdot V_{\mathrm{sub}}\left(t_{\mathrm{0}}\right)+{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}\cdot V_{\mathrm{sub}}\left(t_{\mathrm{6}}\right)\\ & +{\displaystyle \frac{a}{a+b}}\\ & \cdot \mathit{\{}{V}_{\mathrm{std}}\left(t_i\right)-\left({\displaystyle \frac{\mathrm{6}-i}{\mathrm{6}}}\cdot V_{\mathrm{sub}}\left(t_{\mathrm{0}}\right)+{\displaystyle \frac{i}{\mathrm{6}}}\cdot V_{\mathrm{sub}}\left(t_{\mathrm{6}}\right)\right)\mathit{\}}\\ & +{\displaystyle \frac{b}{a+b}}\\ & \cdot \mathit{\{}{V}_{\mathrm{std}}\left(t_j\right)-\left({\displaystyle \frac{\mathrm{42}-j}{\mathrm{6}}}\cdot V_{\mathrm{sub}}\left(t_{\mathrm{36}}\right)+{\displaystyle \frac{j-\mathrm{36}}{\mathrm{6}}}\cdot V_{\mathrm{sub}}\left(t_{\mathrm{42}}\right)\right)\mathit{\}}\end{array}\end{array}$$

$$\begin{array}{}\text{(19)}& \begin{array}{rl}& {\left({\widehat{\mathit{\sigma }}}_{\mathrm{std}\left(t_i,t_j\right)}^{\mathrm{BE}}\left(t_{\mathrm{4}}\right)\right)}^{\mathrm{2}}={\left({\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\cdot {\mathit{\sigma }}_{\mathrm{sub}}\left(t_{\mathrm{0}}\right)\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}\cdot {\mathit{\sigma }}_{\mathrm{sub}}\left(t_{\mathrm{6}}\right)\right)}^{\mathrm{2}}\\ & +{\left({\displaystyle \frac{a}{a+b}}\cdot {\mathit{\sigma }}_{\mathrm{std}}\left(t_i\right)\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{a}{a+b}}\cdot {\displaystyle \frac{\mathrm{6}-i}{\mathrm{6}}}\cdot {\mathit{\sigma }}_{\mathrm{sub}}\left(t_{\mathrm{0}}\right)\right)}^{\mathrm{2}}\\ & +{\left({\displaystyle \frac{a}{a+b}}\cdot {\displaystyle \frac{i}{\mathrm{6}}}\cdot {\mathit{\sigma }}_{\mathrm{sub}}\left(t_{\mathrm{6}}\right)\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{b}{a+b}}\cdot {\mathit{\sigma }}_{\mathrm{std}}\left(t_j\right)\right)}^{\mathrm{2}}\\ & +{\left({\displaystyle \frac{b}{a+b}}\cdot {\displaystyle \frac{\mathrm{42}-j}{\mathrm{6}}}\cdot {\mathit{\sigma }}_{\mathrm{sub}}\left(t_{\mathrm{36}}\right)\right)}^{\mathrm{2}}\\ & +{\left({\displaystyle \frac{b}{a+b}}\cdot {\displaystyle \frac{j-\mathrm{36}}{\mathrm{6}}}\cdot {\mathit{\sigma }}_{\mathrm{sub}}\left(t_{\mathrm{42}}\right)\right)}^{\mathrm{2}}\end{array}\end{array}$$

Figure 3 Schematic diagram of the estimation method for the output of the standard gas at the time of sample air measurement.

2.4.2 Estimation of sample air mole fraction and its combined standard uncertainty using a calibration line

We calculated a calibration line with the estimated outputs of standard gases $\left({\widehat{V}}_{\mathrm{std}\left(t_i,t_j\right)}^{\mathrm{BE}}\left(t_k\right)\right)$ and their variances $\left({\left({\widehat{\mathit{\sigma }}}_{\mathrm{std}\left(t_i,t_j\right)}^{\mathrm{BE}}\left(t_k\right)\right)}^{\mathrm{2}}\right)$ at the time of the sample measurement obtained in Sect. 2.4.1. Although the NDIR output may be regressed with a polynomial equation (e.g., Tanaka et al., 1983), there is no significant difference between the linear and quadratic regression results in this system (shown in Sect. 2.5). To perform a weighted regression with the output of each standard gas, a linear line (y = Sx+I; y: output in voltage, x: mole fraction for CO₂ and log(mole fraction) for CH₄) was adopted for the calibration line (Fig. 4).

Figure 4 Schematic diagram for estimating the CO₂ and CH₄ mole fraction (x) of the gases and the combined standard uncertainty (u_c(x)) from the output in voltage (y) with its SD (σ_y). The gray line indicates the estimated linear calibration line $(y=Sx+I)$.

Following the likelihood method, we identified the slope (S) and intercept (I) for every sample time (k) at the maximum of the likelihood function (L). Solving the normal equation of $\left\{\begin{array}{l}\frac{\partial L}{\partial S}=\mathrm{0}\\ \frac{\partial L}{\partial I}=\mathrm{0}\end{array}\right.$, S and I were obtained as follows.

$$\begin{array}{}\text{(20)}& S\left(k\right)={\displaystyle \frac{\left(\sum w_{ijk}\right)\left(\sum w_{ijk}{x}_i{y}_{ijk}\right)-\left(\sum w_{ijk}{x}_i\right)\left(\sum w_{ijk}{y}_{ijk}\right)}{\left(\sum w_{ijk}\right)\left(\sum w_{ijk}{x}_i^{\mathrm{2}}\right)-{\left(\sum w_{ijk}{x}_i\right)}^{\mathrm{2}}}}\end{array}$$

$$\begin{array}{}\text{(21)}& I\left(k\right)={\displaystyle \frac{\left(\sum w_{ijk}{y}_{ijk}\right)\left(\sum w_{ijk}{x}_i^{\mathrm{2}}\right)-\left(\sum w_{ijk}{x}_i{y}_{ijk}\right)\left(\sum w_{ijk}{x}_i\right)}{\left(\sum w_{ijk}\right)\left(\sum w_{ijk}{x}_i^{\mathrm{2}}\right)-{\left(\sum w_{ijk}{x}_i\right)}^{\mathrm{2}}}}\end{array}$$

Here, x_i is the WS-gas mole fraction determined against the NIES scale. y_ijk is the estimated output of standard gas $\left({\widehat{V}}_{\mathrm{std}\left(t_i,t_j\right)}^{\mathrm{BE}}\left(t_k\right)\right)$ and w_ijk is the reciprocal of the variance $\left(\mathrm{1}/({\widehat{\mathit{\sigma }}}_{\mathrm{std}\left(t_i,t_j\right)}^{\mathrm{BE}}\left(t_k\right){)}^{\mathrm{2}}\right)$. Σ indicates the sum of i (three standard gases) and the same for the following discussion. As shown in Sect. 2.4.1, the combinations of (i,j) are (1,37), (2,38), and (3,39). We omitted i, j, and k for the following expression. The linearity of the calibration line was assessed using the correlation coefficient, and data were rejected when the linearity was deemed insufficient.

The inverse function was used because we estimated the mole fraction from the output in voltage. Furthermore, because the calibration line passes through the weighted mean point $\left(\stackrel{\mathrm{‾}}{x},\stackrel{\mathrm{‾}}{y}\right)=(\frac{\sum wx}{\sum w},\frac{\sum wy}{\sum w})$, practically we used the following line:

$$\begin{array}{}\text{(22)}& x={\displaystyle \frac{y-\stackrel{\mathrm{‾}}{y}}{S}}+\stackrel{\mathrm{‾}}{x}.\end{array}$$

The square of combined standard uncertainty (u_c(x)) for the estimated mole fraction (x) was calculated with the following equation (WMO, 2020):

$$\begin{array}{}\text{(23a)}& \begin{array}{rl}{u}_{\mathrm{c}}^{\mathrm{2}}\left(x\right)& ={\left({\displaystyle \frac{\partial x}{\partial y}}\right)}^{\mathrm{2}}{u}^{\mathrm{2}}\left(y\right)+{\left({\displaystyle \frac{\partial x}{\partial \stackrel{\mathrm{‾}}{y}}}\right)}^{\mathrm{2}}{u}^{\mathrm{2}}\left(\stackrel{\mathrm{‾}}{y}\right)\\ & +{\left({\displaystyle \frac{\partial x}{\partial S}}\right)}^{\mathrm{2}}{u}^{\mathrm{2}}\left(S\right)+{\left({\displaystyle \frac{\partial x}{\partial \stackrel{\mathrm{‾}}{x}}}\right)}^{\mathrm{2}}{u}^{\mathrm{2}}\left(\stackrel{\mathrm{‾}}{x}\right),\end{array}\end{array}$$

where u² is variance for each component. The first term expresses the contribution from the variation in output of the measured air (σ_y) and 60 repeated measurements:

$${\left({\displaystyle \frac{\partial x}{\partial y}}\right)}^{\mathrm{2}}{u}^{\mathrm{2}}\left(y\right)={\displaystyle \frac{{\mathit{\sigma }}_y^{\mathrm{2}}}{S^{\mathrm{2}}}}\cdot {\displaystyle \frac{\mathrm{1}}{\mathrm{60}}}.$$

The second term expresses the contribution from the variation in $\stackrel{\mathrm{‾}}{y}$:

$${\left({\displaystyle \frac{\partial x}{\partial \stackrel{\mathrm{‾}}{y}}}\right)}^{\mathrm{2}}{u}^{\mathrm{2}}\left(\stackrel{\mathrm{‾}}{y}\right)={\displaystyle \frac{\mathrm{1}}{S^{\mathrm{2}}}}{\displaystyle \frac{\sum w^{\mathrm{2}}{\mathit{\sigma }}^{\mathrm{2}}}{{\left(\sum w\right)}^{\mathrm{2}}}}={\displaystyle \frac{\mathrm{1}}{S^{\mathrm{2}}}}\cdot {\displaystyle \frac{\mathrm{1}}{\sum w}},$$

where σ² is the variance of the output for standard gases constituting the calibration line. The third term expresses the contribution from the variation in the slope of the calibration line (S):

$$\begin{array}{rl}& {\left({\displaystyle \frac{\partial x}{\partial S}}\right)}^{\mathrm{2}}{u}^{\mathrm{2}}\left(S\right)\\ & ={\displaystyle \frac{{\left(y-\frac{\sum wy}{\sum w}\right)}^{\mathrm{2}}}{S^{\mathrm{4}}}}\cdot \sum {\mathit{\sigma }}^{\mathrm{2}}{\left({\displaystyle \frac{\partial S}{\partial y}}\right)}^{\mathrm{2}}\\ & ={\displaystyle \frac{(y-\frac{\sum wy}{\sum w}{)}^{\mathrm{2}}}{S^{\mathrm{4}}}}\cdot {\displaystyle \frac{\sum w}{\left(\sum w\right)\left(\sum w{x}^{\mathrm{2}}\right)-{\left(\sum wx\right)}^{\mathrm{2}}}}.\end{array}$$

The fourth term expresses the contribution from the variation in $\stackrel{\mathrm{‾}}{x}$. The NIES 09 CO₂ scale is based on the gravimetric primary standard gases using a one-step dilution (Machida et al., 2011). When 15 gravimetrically prepared mixed gases in the range of 350–390 ppm were measured on an existing scale, the SD of the differences from their mole fractions determined on that existing scale was 0.042 ppm (Tohjima et al., 2006), which we adopt as the uncertainty of the gravimetric method. For the NIES 09 CO₂ scale, we prepared eight standard gases gravimetrically in the range of 340–450 ppm for atmospheric measurements. However, we maintain as primary standards of the NIES 09 CO₂ scale the calibration results of eight different standard gases in the same concentration range that had demonstrated long-term stability, calibrated using the gravimetrically prepared standards (Machida et al., 2011). During calibration of primary standards, we performed repeated measurements (N = 27–40) and used their mean values as calibration values. The SD of these measurements was 0.01–0.02 ppm (similar values were observed for calibrations of secondary and working standards). Since the mole fractions of standards prepared by the one-step gravimetric method are independent, the uncertainty for each primary standard becomes 0.02 ppm $\left(\sqrt{\frac{{\mathrm{0.042}}^{\mathrm{2}}}{\mathrm{8}}+{\text{SD}}^{\mathrm{2}}/N}\right)$, including the transfer uncertainty. In effect, as the propagation term is an order of magnitude smaller, the root of the fourth term in Eq. (23a) also becomes 0.02 ppm.

On the one hand, the NIES 94 CH₄ scale is based on the gravimetric primary standard gases using a four-step dilution. The root of the fourth term is estimated to be 3.2 ppb (the Supplement). These values of the fourth term are common to all data points, so they should be denoted separately. Summarizing the first three terms, the u_c(x) for the estimated mole fraction (x) is as follows.

$$\begin{array}{}\text{(23b)}& \begin{array}{rl}& u_{\mathrm{c}}\left(x\right)={\displaystyle \frac{\mathrm{1}}{S}}\\ & \cdot \sqrt{{\displaystyle \frac{{\mathit{\sigma }}_y^{\mathrm{2}}}{\mathrm{60}}}+{\displaystyle \frac{\mathrm{1}}{\sum w}}+{\displaystyle \frac{\mathrm{1}}{S^{\mathrm{2}}}}\cdot {\displaystyle \frac{(\left(\sum w\right)y-\sum wy{)}^{\mathrm{2}}}{{\left(\sum w\right)}^{\mathrm{2}}\left(\sum w{x}^{\mathrm{2}}\right)-\left(\sum w\right){\left(\sum wx\right)}^{\mathrm{2}}}}}\end{array}\end{array}$$

Since the calculation was done for the logarithm of the mole fraction for CH₄, the u_c(x) for the estimated mole fraction was determined differently for the higher level $(u^{+}=x\left(e^{u_{\mathrm{c}}\left(x\right)}-\mathrm{1}\right))$ and lower level $(u^{-}=x\left(\mathrm{1}-e^{-u_{\mathrm{c}}\left(x\right)}\right))$. However, the average value is expressed as the u_c(x) since the difference is less than 0.1 ppb in real terms.

Figures S3 to S11 show the time series of the u_c(x) for the ambient air CO₂ mole fraction. Most uncertainties were distributed around 0.05 ppm, but they can be higher than 0.3 ppm, especially during summer. Note that the SD of the output (σ_y) of the sample air could become significant due to large diurnal variation during summer (Sasakawa et al., 2013) since the output of the sample air could include a natural variation of the atmosphere during the measurement period of 3 min. Therefore, following the recommendations outlined in the “Sensor precision and atmospheric variability” chapter of the WMO (2020), it is valuable to concurrently present the sensor repeatability (u_r). The sensor repeatability should ideally be determined by conducting continuous measurements with a cylinder connected to the sample inlet. However, this approach has been impractical since the system's initial installation. Nevertheless, in our measurement system, SWS-gas measurements are conducted before and after each sample measurement, whereby constant gas from the cylinder is measured continuously for 3 min. We estimated u_r by calculating the u_c(x) for the SWS-gas following procedures analogous to those used for sample calculations and then temporally interpolating these values between the pre- and post-sample SWS-gas measurements over the sample measurement period. The estimated u_r values were distributed so as to closely follow the minimum values of sample u_c(x) (Figs. S3–S11). In rare cases, they were slightly higher than the sample minimum values, which occurred when the mole fraction of the SWS-gas exceeded that of the highest WS-gas, resulting in increased uncertainty from the calibration curve.

Figures S12 to S17 show the time series of the u_c(x) for the ambient air CH₄ mole fraction. Most are within 5 ppb, but they can be above 10 ppb during summer at KRS and DEM. Meanwhile, at NOY, although the baseline is below 5 ppb, values exceeding 10 ppb can occur regardless of season. All sites are affected by short-term variations from summer wetland emissions, while NOY is presumed to experience additional significant short-term variations from anthropogenic methane sources such as leakage of natural gas. The u_r was estimated using a method similar to that for CO₂, and the values were distributed so as to closely follow the minimum values of sample u_c(x). Although values remained below 5 ppb, they could fluctuate from near 0 to 5 ppb over several days. Such fluctuations were also apparent in the samples and are considered to be mainly due to sensor stability.

2.5 Stability check with the SWS-gas measurement

We calculated a calibration line only when the SWS-gas measurements closest to both sides of the sample measurements were normal. The normality of the SWS-gas measurement was assessed as follows. The same on-site compressed air was measured several times (for about a week) since the air was used as an SWS-gas until its pressure dropped below 0.1 MPa. The on-site compressed air output value would vary smoothly if the system were stable. When the system temporarily became unstable the corresponding large changes of the analyzer output could not be corrected by the SWS to a sufficient degree. To identify such occasions, we first estimated the output of standard gases at the time of the target SWS-gas measurement by interpolating Δ^B(t_i) and Δ^E(t_j) based on the output value of the target SWS-gas itself (Sect. 2.4.1). Then, calculating a calibration line with the estimated output of standard gases, we obtained the mole fraction of the target SWS-gas. Second, we estimated the output value of the SWS-gas at the time of the target SWS-gas measurement by interpolating the outputs of the two adjacent SWS-gases. Then, we determined the mole fraction of the target SWS-gas in the same manner. If these estimated mole fractions differed from each other by more than 1 ppm for CO₂ and 10 ppb for CH₄, we regarded the target SWS-gas data as abnormal. This assessment (referred to as a “self-check-value (scv) for SWS-gases” in Fig. 5) was done while the adjacent SWS-gas measurements were conducted for the same on-site compressed air.

We then checked the system's stability with the measurement of the SWS-gas. Interpolating the outputs of the SWS-gases adjacent to the standard gas measurements, we calculated the mole fractions of the SWS-gas with the calibration line at the time of STD2, which was used to assess the coefficient of determination in Sect. 2.3. We regarded the estimated mole fractions of the SWS-gas as independent; thus, we obtained 14 estimated mole fractions if the measurements for the same SWS-gas continued for a week. We determined a threshold for the SD (σ_sws; 1 ppm for CO₂, 10 ppb for CH₄) and the fluctuation range (3 ppm for CO₂ and 30 ppb for CH₄) obtained from the estimated independent dataset. All the data that exceeded the threshold were deleted. A flowchart for the calculation method is shown in Fig. 5.

Figure 5 Flowchart of the mole fraction calculation method.

3 Reproducibility

3.1 Change in temporal conditions

The σ_sws values that remained after meeting the criteria described in Sect. 2.5 serve as an indicator of measurement reproducibility at 12 h intervals throughout the given period (often a week). Regardless of the site or time of measurement, σ_sws for CO₂ mostly was below 0.2 ppm, and it was below 5 ppb for CH₄ (Figs. S3–S17). After introducing CRDS, the working standard gas measurement interval was changed to 48 h. At the same time, since the consumption flow rate doubled, the SWS-gas changeover time decreased to approximately 3 d, which makes the number of SWS measurements with NDIR/TOS only three at most. Although the σ_sws values obtained from a few data points are for reference only, they were distributed in almost the same range (light blue dots in Figs. S3–S17).

3.2 Change in sensor and temporal conditions

The CRDS instrument operated at KRS from July 2015, DEM from June 2016, and NOY from August 2016, albeit for a short period of time within JR-STATION’s long observation period. The CRDS instrument is a highly stable analyzer used in greenhouse gas observations worldwide (Yver Kwok et al., 2015). We compared the recalculated NDIR (and TOS) values with the CRDS values to check the long-term reproducibility (Figs. 6–7, S18–S21). The flow path branches off after the six-port valve (Fig. 1), so the same air is analyzed. The CRDS instrument operates independently of the existing system, so information on instrument error flags and valve switching timing is not linked to the measured data. Therefore, the timing of the standard gas measurement was captured by detecting the CH₄ mole fraction of the lowest standard gas. Only data from periods were extracted when the SWS results fulfilled the criteria outlined in the previous section. The CRDS output values were converted to the NIES scale based on the WS measurements and averaged over 3 min for comparison. The temperature in the warm box of the CRDS instrument (data column name is “WarmBoxTemperature”) was kept constant (45.00 °C), but it may vary significantly in rare cases. CRDS outputs data every few seconds, but the output frequency may drop abnormally. Since the device is considered unstable under these conditions, the data were not used for comparison if the temperature change was more than 0.03 °C or if the number of data points was fewer than 10 in 3 min averaged over the CRDS data. Since some observed values fell outside the mole fraction range of the standard gas, only values within the range were used to calculate the difference between the two.

Figure 6 Relationship between the CO₂ mole fraction by CRDS at KRS and the difference in the respective CO₂ mole fractions measured by the NDIR analyzer and CRDS (NDIR − CRDS): the CRDS values were averaged over the corresponding 3 min period. The light blue (gray) circle shows the difference from a low-altitude (high-altitude) inlet. The dotted lines indicate the mole fraction of standard gases. The panel (mean ± SD) in the bottom right represents the average difference for each inlet. Only data that were within the standard gas mole fraction range were used. The red dots indicate the values averaged every 1 ppm for the combined high-altitude and low-altitude data. However, calculations were only made when the number of data points used was 100 or more. The results from the calibration curve using the quadratic equation are shown as green dots. Error bars indicate the SD. These annual averages are marked by their respective colors in the upper left.

Figure 7 Relationship between the CH₄ mole fraction by CRDS at KRS and the difference in the respective CH₄ mole fractions measured by the TOS and CRDS (TOS − CRDS): the CRDS values were averaged over the corresponding 3 min period. The light blue (gray) circle shows the difference from a low-altitude (high-altitude) inlet. The dotted lines indicate the mole fraction of standard gases. The panel (mean ± SD) in the left bottom represents the average difference for each inlet. Only data that were within the standard gas mole fraction range were used. The red dots indicate the values averaged every 10 ppb for the combined high-altitude and low-altitude data. However, calculations were only made when the number of data points used was 100 or more. Error bars indicate the SD. The annual average is marked in the upper right.

As for the NDIR analyzer, there was no significant difference between the high-inlet and low-inlet differences, indicating no bias due to differences in inlets (Figs. 6, S18–S20). However, regardless of the year of observation, the NDIR analyzer showed lower values than those of CRDS by about 0.1 ppm at all sites. The CO₂ mole fraction has more significant diurnal variation during the summer months (Sasakawa et al., 2013), so the error bars were more prominent (Fig. S20). Still, the amount of bias remains the same. As mentioned in Sect. 2.1, using our isotopically lighter standard gas, the NDIR analyzer measures a lower CO₂ mole fraction in sample air with an actual CO₂ mole fraction. This 0.1 ppm difference can be attributed to the optical filter property of the NDIR analyzer used in this system. In addition, this bias does not change over time, indicating that this system produces results with good reproducibility over time. We plan to make the correction of 0.1 ppm for the NDIR analyzer in the published data.

Results from data using calibration curves with quadratic equations are also shown. There were no significant differences from the results using the linear regression in any concentration range, confirming the linearity of the NDIR analyzer used in this system. In the first place, the coefficients of the second order of the regressed quadratic equation frequently changed, both positively and negatively, so it is unlikely that the essential response of NDIR can be estimated by a polynomial equation in the concentration range of the standard gas used in this system.

As for the TOS, there was also no significant difference between the high-inlet and low-inlet differences, indicating no bias due to differences in inlets (Fig. 7). As an overall average within the actual measurement range in Siberia, the TOS did not differ from CRDS, but there can be a bias from −10 ± 5 to 5 ± 5 ppb with the CRDS output, depending on the mole fraction. However, the degree of the bias varied from year to year, and its cause was unknown, so a constant correction cannot be made. The time series showed higher values during the winter period in some years (Fig. S21), which may be the result of the increase in the mole-fraction-dependent difference seen in Fig. 7.

At NOY and DEM, it was discovered that the temperature controller of the catalytic unit was not functioning correctly. Since the TOS is sensitive to CO and H₂ in the air, it could produce unusually high values without a proper catalytic unit. For the period, only the data from CRDS should be published. Since no catalytic unit errors were identified at the other sites, the ambient atmospheric values were detected, as is the case with KRS.

Both the NDIR analyzer and the TOS have shown good reproducibility over several years for CRDS, suggesting that this measurement is convincing.

4 Full uncertainty estimates

In Table 3, we present all uncertainties estimated in this study. Since our measurements are based on gravimetrically prepared primary standards, we have also calculated the uncertainties propagated from the preparation of primary standards through to working standard gases (see Sect. 2.4.2 and the Supplement). Our observation system records 3 min average mole fractions of atmospheric CO₂ and CH₄ at multiple heights every hour. These average values are published with their u_c(x). While our system corrects for sensor drift using on-site air (SWS-gas), we estimated the sensor repeatability (u_r) during 3 min sample measurements using the u_c(x) of 3 min measurements of this SWS-gas (Sect. 2.4.2). The SWS-gas was used for approximately 1 week, and its reproducibility was estimated from the repeated measurements (Sect. 3.1). We verified reproducibility for three sites where CRDS was installed by comparing the CRDS and NDIR (or TOS) measurements (Sect. 3.2). For CO₂, the NDIR analyzer showed 0.1 ppm lower values, which we attributed to isotope effects in NDIR; consequently, we plan to implement corrections for NDIR measurements. For CH₄, annual mean values agreed within the measurement scatter. However, within the range of mole fractions measured each year, differences from −10 ± 5 to 5 ± 5 ppb occurred.

Table 3 Summary of uncertainties.

^* Isotope effect in the NDIR analyzer. We plan to apply this correction to the NDIR measurements.

5 Conclusions

Using the standard-gas-saving system, we have developed and validated an improved methodology for calculating mole fractions and uncertainties in long-term CO₂ and CH₄ measurements at remote sites. Following the recommendations shown in the GAW report (WMO, 2020), we provided both the combined standard uncertainty (u_c(x)) and sensor repeatability (u_r) for each measurement, enabling data users to select appropriate data based on their specific research requirements. Comparison with CRDS measurements demonstrated excellent long-term reproducibility of the NDIR analyzer and TOS, with the NDIR analyzer showing a consistent offset due to isotope effects that can be corrected. For CH₄, while annual means agreed well with CRDS, some differences were observed depending on the year and mole fraction level. The working standard-gas-saving system, using on-site compressed air as a sub-working standard gas, proved effective for long-term monitoring at remote sites, showing good temporal reproducibility. By quantifying all major uncertainty components, including those from primary standard gas preparation to working standards, we have provided a comprehensive uncertainty budget that enhances the value of the JR-STATION dataset for studying greenhouse gas dynamics in Siberia. The validated performance of our cost-effective system design demonstrates its suitability for maintaining high-quality, long-term greenhouse gas monitoring networks in remote regions where infrastructure is limited.