A new procedure for processing eddy-covariance data to better 1 quantify atmosphere-aquatic ecosystem CO 2 exchanges

1 quantify atmosphere-aquatic ecosystem CO2 exchanges 2 3 Tatsuki Tokoro1,2 Tomohiro Kuwae1 4 1Coastal and Estuarine Research Group, Port and Airport Research Institute, 3-1-1 5 Nagase, Yokosuka 239-0826, Japan 6 2Present address: Marine Macrophyte Ecosystem Group, National Research Institute of 7 Fisheries and Environment of Inland Sea, 2-17-5, Hatsukaichi 739-0452, Japan 8 Correspondence to: T. Tokoro (tokoro-t@ipc.pari.go.jp) 9 10 Abstract. The capture of carbon by aquatic ecosystems and its sequestration in 11 sediments has been studied as a potential method for mitigating the adverse effects of 12 climate change. However, the evaluation of in situ atmospheric CO2 fluxes is 13 challenging because of the difficulty in making continuous measurements over areas 14 and for periods of time that are environmentally relevant. The eddy covariance (EC) 15 method is the most promising approach to address this concern with the measurement of 16 atmospheric CO2 fluxes. However, methods to process the data obtained from EC 17 measurements are still being developed, and the estimated air-water CO2 fluxes have 18 large uncertainties and differ from those obtained using conventional methods. In this 19 study, we improved the post-processing procedure for the EC method to reduce the 20 uncertainty in the measured air-water CO2 fluxes. Our new procedure efficiently 21 removes erroneous fluxes using a combination of filtering methods based on the 22 received signal strength indicator of the EC sensor, the normalized standard deviation of 23


Introduction
Aquatic environments are considered critical to the mitigation of adverse climate change effects because of their ability to store atmospheric CO2.Previous studies have estimated that the ocean absorbs approximately one-fourth of the CO2 emitted by anthropogenic activities (IPCC, 2013).However, the effect of shallow aquatic ecosystems on atmospheric CO2 remains a controversial topic.Several previous studies have concluded that shallow aquatic ecosystems are sources of atmospheric CO2 after taking account of carbon inputs from land (e.g.Gazeau et al., 2005;Borges et al., 2006;Chen et al., 2013).In contrast, some autotrophic, shallow aquatic ecosystems have been reported to be net sinks for atmospheric CO2 (e.g.Schindler et al., 1997;Tokoro et al., 2014).
In situ measurements of atmospheric CO2 fluxes are necessary for precise analysis Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-499Manuscript under review for journal Biogeosciences Discussion started: 14 December 2017 c Author(s) 2017.CC BY 4.0 License.Merlivat, 1986;Wanninkhof, 1992;Ho et al., 2006).In the case of shallow systems, water velocity fields and depths have been used to estimate the gas transfer velocity (O'Conner and Dobbins, 1958;Borges et al., 2004).In this study, the method using such empirical model for estimating aquatic CO2 flux is defined as the Bulk Formula method (BF method).
However, application of the BF method is limited because of its poor temporal and spatial coverage.Moreover, in previous studies, air-water CO2 fluxes have been determined mostly as snapshots that did not account for diurnal changes or annual cycles, the result being considerable uncertainty and bias (Kuwae et al., 2016).In brackish environments in particular, temporal variability of water fCO2 is significant, whereas the carbonate buffer effect is weak, and the fluctuations of fCO2 become very large (Zeebe and Wolf-Gladrow, 2001).Use of BF methods to carry out a comprehensive analysis of dynamic carbon cycling in aquatic environments with large spatial and temporal variability would therefore be very costly and require much effort.
Another method for evaluating the air-water CO2 fluxes is direct measurement of in situ fluxes, which involves use of a chamber floating on the water surface (e.g.Frankignoulle, 1988;Tokoro et al., 2008) and eddy covariance devices (vide infra).The floating chamber method is used to determine the air-water CO2 flux from continuous measurements of CO2 concentrations in the air inside a hollow, box-shaped device floating on the water surface.Although this method is the easiest of the direct methods to use in shallow coastal waters because of its relative simplicity, like the BF method, it is poorly suited for obtaining long-term measurements over wide areas.
The eddy covariance (EC) method, which is commonly used to determine mass and heat fluxes in terrestrial environments, has recently been used to estimate air-water fluxes of greenhouse gases (e.g. Lee et al., 2004).The determination of the EC CO2 flux is based on the micro-meteorological behavior of atmospheric eddy diffusion and is calculated from the covariance of atmospheric CO2 concentrations and vertical wind speeds measured at high frequency (more than 10 Hz).Because EC measurements can be performed automatically and represent the flux over a large area, the EC method can be used to obtain a detailed analysis of CO2 fluxes.
Despite the promise of EC measurements, application of the EC method in aquatic environments remains challenging (Tsukamoto et al., 2004;Rutgersson and Smedman, 2010;Vesala, 2012;Blomquist et al., 2013;Kondo et al., 2014;Ikawa and Oechel, 2014;Landwehr et al., 2014).The difficulty of making aquatic EC measurements is that the air-water CO2 flux is small compared with the air-land CO2 flux (Vesala, 2012;Landwhehr et al., 2014).The main technical problem with the EC method is cross sensitivity, which reflects the interference between the atmospheric CO2 and H2O measurements caused by spectrometric error (Kohsiek et al., 2000;Prytherch et al., 2010;Kondo et al., 2014;Landwehr et al., 2014).A procedure based on the relationship between atmospheric CO2 concentration and relative humidity called the PKT correction has been proposed to correct for the effects of cross sensitivity (Prytherch et al., 2010).However, the PKT correction is not always effective.Past studies have revealed that cross sensitivity is reduced only after application of certain operational procedures such as cleaning the optical lens of the sensor (Ikawa and Oechel, 2014;Kondo et al., 2014) and drying the sample gas (Landwhehr et al., 2014).
There are several problems in addition to cross sensitivity in using EC measurements in aquatic environments.The uncertainty of EC measurements has been attributed to the spatial and temporal heterogeneity of water (Mørk et al., 2014).The EC flux is calculated as the average within a measurement area called the "footprint", which can range from several hundred meters to several kilometers windward from the measurement point (Schuepp et al., 1990).Therefore, EC fluxes at heterogeneous water sites are different from the fluxes determined by methods that estimate the CO2 flux in an area of only several square meters (e.g., the BF method and floating chamber method).The EC flux is an average flux over a certain time interval (approximately several tens of minutes) (Lee et al., 2004), whereas the BF method estimates the flux at the time of sampling.Thus, EC fluxes estimated at sites where fluxes are temporally variable also differ from fluxes obtained using other methods.Furthermore, the inflow of terrestrial air into the measurement site can generate uncertainty in the flux measurement because the atmospheric CO2 concentration over terrestrial vegetation may differ significantly from the concentration over water.The inflow of terrestrial air can cause unnatural temporal changes in the atmospheric CO2 concentration and spatial heterogeneity at the measurement site.It is therefore necessary to account for the characteristics of the aquatic environment and apply appropriate post-processing (PP) procedures (Leinweber et al., 2009) to avoid large uncertainties or biases in EC flux calculations.
In this study, we developed a PP procedure for EC aquatic measurements.This PP procedure involves the exclusion of erroneous data and correction of unnatural changes in the atmospheric CO2 using a series of data-filtering steps.The new process is based on the idea that cross sensitivity and environmental heterogeneity during flux measurement cause spikes, drifts, offsets, and long-term variation in the CO2 and H2O raw data.We compared the results calculated with our new PP procedure to those obtained using conventional EC PP procedures along with BF flux data as an example Page 6 of 40 Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-499Manuscript under review for journal Biogeosciences Discussion started: 14 December 2017 c Author(s) 2017.CC BY 4.0 License. of a frequently used method.We then discuss the differences between the conventional BF and EC fluxes with respect to atmospheric and environmental characteristics.

Field measurements
Continuous EC measurement data were used for the evaluation of the PP procedure and the analysis of atmospheric-aquatic ecosystem CO2 exchange.The data were collected from a brackish lagoon in Japan (the Furen Lagoon, Fig. 1) from 28 May to 21 October 2014, during which time the water surface was not frozen.Most of the study area (57.4 km 2 ) is covered by seagrass meadows (mainly Zostera marina).The water is shallow (1-2 m), except in a channel that connects the eastern and western basins of the lagoon (depth = approximately 5 m).Freshwater flows into the western basin through several rivers that run through the surrounding grass farms, and seawater is exchanged through the lagoon mouth, which opens to the Okhotsk Sea.A previous study has found that the air-water CO2 flux in the lagoon is affected by changes of salinity caused by the inflow of river water and tides as well as by changes of dissolved inorganic carbon resulting from biological processes such as photosynthesis (Tokoro et al., 2014).The measurement platform was built at the same site used in the study of Tokoro et al. (2014) (N43° 19.775', E145° 15.463'); the effects of photosynthesis and changes in salinity are most notable at this location in the lagoon (Tokoro et al., 2014).
The EC devices used in this study were as follows.Atmospheric CO2 concentrations and water vapor were measured with an open-path sensor (LI-7500A, LI-COR, USA).
The three-dimensional (3D) wind velocity, air temperature, and atmospheric pressure were measured with a 3D acoustic Doppler velocimeter (CSAT-3, Campbell Scientific, Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-499Manuscript under review for journal Biogeosciences Discussion started: 14 December 2017 c Author(s) 2017.CC BY 4.0 License.USA).The data were logged and managed by a SMARTFlux system (LI-COR).The open-path sensor and the wind velocimeter were attached to the platform approximately 3.0-5.5 m above the water surface (the height varied with the tide).The sampling rate for all data was 10 Hz, and the fluxes (CO2, water vapor, and heat) were calculated as averages over 30-min intervals.Batteries and solar panels were attached to the platform as a power source.Battery replacement, data collection, and device maintenance were performed approximately every two weeks.Water temperature and salinity were measured continuously with a conductivity-temperature sensor (Compact-CT, Alec, Japan).

Calculation of fluxes using the conventional PP procedure (PP1)
The air-water CO2 flux (F) was calculated every 30-min using the following equation: , where the coefficients F1 and F2 are correction terms based on the transfer functions that correct for the frequency attenuation of the air-sea CO2 flux caused by the response time of the sensor, path-length averaging, sensor separation, signal processing, and fluxaveraging time (Massman, 2000).The first term on the right-hand side of Eq. ( 1) is the product of F1 and the uncorrected air-sea CO2 flux calculated as the covariance of the CO2 density ρc and the vertical wind speed w (the bar and the prime indicate the mean and the deviation from the mean, respectively).The second and third terms are the Webb-Pearman-Leuning correction of latent heat and sensible heat, respectively (Webb et al., 1980).The other variables in Eq. ( 1) are defined as follows: ρd = dry air density; ρv = water vapor density; Ta = air temperature; and μ = ratio of the molar weight of dry air to that of water vapor.The wind speed was corrected by using a double rotation to make the average vertical wind speed zero during the 30-min time interval (Lee et al., 2004).The footprint (measurement area) depends on several factors, including the measurement height, wind speed, atmospheric stability, and measurement site roughness (10 −4 cm) (Schuepp et al., 1990).This footprint ranged from several hundred square meters to several square kilometers on the windward side of the measurement site.
The deviation of each parameter in Eq. ( 1) from the 30-min average was calculated by subtracting the 30-min average from the instantaneous data after deleting obviously erroneous data (e.g., negative values of CO2 or water vapor concentration).Other corrections to the raw data included coordinate rotation of the 3D wind component (double rotation; Lee et al., 2004), time lag of the measurement due to the separation of the CO2 sensor and the wind velocimeter (covariance maximization; Lee et al., 2000), exclusion of wind data contaminated by the wind velocimeter flame, and correction of the measurement noise (statistical tests of Vickers and Mahrt, 1997) based on the default settings of the data management software (EddyPro 5.1.1,LI-COR).

Calculation of flux using our new PP procedure (PP2)
After calculating the EC flux using the conventional PP procedure (PP1) described in Sect.2.2, we recalculated the EC flux using our new procedure (PP2; Fig. 2).The PP2 procedure is mainly based on filtering and excluding erroneous data rather than on the PKT method of data correction.The PP2 procedure is also focused on aquatic environments in which the spatial and temporal variations in atmospheric CO2 are large.The RSSI is one of the data concurrently obtained from the CO2 sensor of the EC measurement instrumentation and indicates the available signal strength of the sensor.This parameter has been used for evaluating measurement validation of objective gases such as methane.In this study, we used the RSSI to filter the CO2 data because both CO2 and methane absorb infrared radiation.First, data in the 30-min time series were excluded if their RSSI was low.The RSSI threshold for exclusion was set to 90 %, because this percentage was sufficiently high to identify valid CO2 concentrations and allowed most of the measurement data (~81 %) to be retained.This 90 % threshold was also recommended as the quality criterion for the LI-7700 methane analyzer (LI-COR), which has been used to measure methane concentrations.Finally, HP filtering was applied to the calculation of the deviations from the mean of each parameters in Eq. ( 1).This procedure corrected relatively long-term (several minutes to 30 minutes) variations in CO2 or water vapor concentrations that were independent of eddy fluctuations and were caused by the temporal and spatial heterogeneity of the atmospheric mass.HP filtering is often applied to measurements in a complex environment; however, incorrect application of HP filtering results in underestimation of fluxes (Lee et al., 2004).HP filtering was applied by using an exponential moving average as follows: where xi and xi' are an instantaneous datum and deviation from the mean at time i, respectively.The parameter ɤ is the time constant of the exponential moving average, which was determined to be 150 s in a previous study (McMillan, 1988), and f is the sampling frequency (10 Hz).HP filtering was applied to all of the measured data (i.e., 3D wind velocity, air temperature, CO2 and water vapor concentrations, and atmospheric pressure).

Measurement of fluxes using the BF method
BF flux measurements were performed during the daytime on 29 May, 15 July, and 21 September 2014 for comparison with the EC measurements.In the BF method, flux is calculated as follows: where k is the transfer velocity, which was calculated using one or another of the following three empirical equations.The first was the equation of Wanninkhof (1992) (W92).Because this equation was constructed using tracer methodology under oceanic conditions, it might be inappropriate to our measurements because of the differences in fetch and water depth.However, this equation is very commonly used in a variety of oceanographic CO2 flux studies, including coastal measurements.We therefore used it for comparison.The second equation was that of Borges et al. (2004) (B04), which has been applied in an estuarine study that involved use of the floating chamber method.This equation makes use of the current velocity in addition to the wind speed.The third equation was that of Mørk et al., (2014) (M14), which was formulated to characterize the transfer velocity in a fjord for use with the EC method.Thus, the second and third equations can be used to characterize coastal gas transfer velocity.The wind speed for the gas transfer velocity was measured by the EC device and normalized to a height of 10 m from the water surface using a logarithm law (Kondo, 2000).The current speed and depth were measured by the sensors attached to the platform (Compact-TC and Compact-EM, Alec, Japan).The final CO2 flux was calculated as the weighted average of the fluxes measured within the EC footprint (Schuepp et al., 1990).
The parameter S is the dissolution coefficient of CO2, which was estimated from the water temperature and salinity (Weiss, 1974).The parameters fCO2water and fCO2air are the fugacity of water and atmospheric CO2, respectively.Water temperature and salinity were measured with a handheld conductivity-temperature-depth sensor (ACTD-DF, JFE Advantech, Japan).The water samples used to determine fCO2water were collected just below the water surface (up to 20 cm below the water surface) to measure the concentration of CO2 where direct gas exchange with air occurs.The sampling was performed within the EC footprint to minimize the effect of spatial heterogeneity when comparing the BF and EC fluxes.The sampling points were determined from the wind direction and the distance from the platform measured using a hand-held GPS unit (Venture HC, Garmin, USA; see Table S1).The water fCO2 was determined from the total alkalinity and the dissolved inorganic carbon content of the water sample using a batch-type carbonate measurement system (ATT-05, Kimoto electrics, Japan) and the CO2SYS program (Pierrot et al., 2006).

PP1 data
During the measurement period, 4464 flux data points corresponding to 2232 hours were obtained; 1971 of those data points (44 %) were excluded as erroneous data after PP1 application.The mean and SD of the EC CO2 fluxes were −1.93 and 52.4 μmol m −2 s −1 , respectively.Figure 3 Spikes and discontinuities were observed in the atmospheric CO2 and water vapor concentrations, despite the prior de-spiking process applied by the data management system.The cumulative covariance indicated that the covariance at certain periods (0-5 min) contributed significantly to the total cumulative covariance.

PP2 data
The EC CO2 flux data subjected to PP2 (RSSI, nSD, and HP filtering) are shown in Fig. absolute values of some fluxes.Figure 7 shows an example of the results of filtering atmospheric CO2 concentrations as well as the cumulative covariance of the atmospheric CO2 concentrations and vertical wind speed (measured at 8:00 on August 21 [day 84]).These data were not excluded by the RSSI and nSD filtering (RSSI = 100 %, nSD of CO2 = 9.84 × 10 −3 , nSD of H2O = 2.07 × 10 −2 ).The concentration of atmospheric CO2 showed a trend over the 30-min time interval, the indication being that the block average could not extract appropriate eddy movements from the time-series data.The cumulative covariance of the block averaging (BA) was unusually large; after HP filtering, the values were more reasonable because filtering successfully excluded deviations caused by the variation for about 10-min of atmospheric CO2 concentrations.

BF data
The measured BF fluxes showed spatial and seasonal variations (see Table S1).The means and SDs obtained by using the three gas transfer velocity equations to analyze the weighted averages of the BF CO2 flux data were 0.44 ± 0.33 μmol m −2 s −1 (n = 15) on 29 May (day 1), 2.10 ± 1.58 μmol m −2 s −1 (n = 14) on 15 July (day 48), and −0.11 ± 0.13 μmol m −2 s −1 (n = 9) on 21 September (day 115).The differences in the number of data reflect problems during the measurements, such as inadequate water depth and problems with the water-depth sensor.Except for the measurements on 15 July, the difference between the results obtained with the three gas transfer velocity equations were not significant.
The exclusion of erroneous outliers by PP2 also contributed to the time-series analysis.The power spectra of the EC CO2 fluxes after PP1 evidenced large, noise-like fluctuations at high frequencies (Fig. 8), and thus any suggestion of peaks in the time series was obscured.After PP2, however, the noise-like fluctuations were smaller, and two peaks associated with semi-diurnal (~12.5 h) and diurnal (~24 h) time intervals were apparent.The fCO2 variations in the lagoon, which are among the parameters that regulate air-water CO2 fluxes, have been confirmed to be related to mixing of lagoon water with freshwater coming from rivers and with biological processes such as photosynthesis (Tokoro et al., 2014).Given that the former and latter phenomenon are caused by the semi-diurnal tidal cycle and diel changes of irradiance, respectively, the peaks in the power spectra are consistent with the results of Tokoro et al. (2014).This consistency is a good demonstration of the utility of the PP2.
The EC data that were much different from the BF results were excluded by the PP2 (Fig. 9).The remaining EC fluxes in May and September seemed to agree well with the BF fluxes.We believe that one of the reasons that the results were comparable was the improvement in the accuracy of the EC fluxes with the use of PP2.Another reason was our strategy of making BF measurements at multiple points within the EC footprint to This similarity indicates that the effect of currents in the B04 and the effect of earlybreaking waves in the M14 were overestimated at our measurement site.
However, the EC fluxes estimated with PP2 did not always agree with the BF fluxes.Because the fCO2water is theoretically never negative, a theoretical maximum negative BF flux can be calculated by arbitrarily setting fCO2 equal to zero.The maximum negative flux calculated in this way (with M14) was −6.16 μmol m −2 s −1 at 15:00 on 30 May (day 2), when the maximum wind speed was recorded (11.9 m/s).
Forty-seven EC flux data points (3 % of all data) indicated even lower fluxes.
Moreover, the mean of the EC fluxes with PP2 (-0.54 μmol m −2 s −1 ) was more negative or almost the same as the mean of the theoretical BF fluxes (-0.26 μmol m −2 s −1 in, -0.55 μmol m −2 s −1 , and -0.43 μmol m −2 s −1 estimated with W92, B04, and M14, respectively).Because actual BF fluxes include less negative as well as positive fluxes, these EC fluxes cannot be explained by only the BF fluxes.
Similar inconsistencies between air-water CO2 fluxes calculated with the EC method and other conventional methods have been reported in several studies (e.g., Tsukamoto et al., 2004;Rutgersson et al., 2010).In the case of coastal measurements, water-side convection due to the vertical temperature difference inside water has been postulated to enhance the gas transfer velocity (Rutgersson et al., 2010).However, such an enhancement was not previously observed with direct flux measurement using a floating chamber at our site (Tokoro et al., 2014).Because of the very shallow water depth (less than 2 m) at our site, we suspect that water-side convection was weak and was not the main reason for the inconsistency of the fluxes.Assuming that the EC fluxes obtained with PP2 are valid, the discrepancy between the EC and BF fluxes was also postulated to reflect the limitations of the BF method and/or the difference of the measurement height of the BF and EC methods; the former and the latter are at the water surface and the height of the EC devices, respectively.As for the BF method limitaitons, seagrass leaves, which reached the water surface during low tide at the study site, might have affected the physical and chemical conditions at the water surface (Watanabe and Kuwae, 2015).In BF theory, the CO2 flux is caused by the CO2 concentration gradient just below the water surface.The BF method should therefore not be applied when seagrass is present on the water surface.A previous study that investigated the radiocarbon isotopic signatures of seagrass at the study site indicated that of the total CO2 assimilated by the seagrass, 0-40 % (mean = 17 %) originated from the atmosphere and the rest from the water (Watanabe and Kuwae, 2015).The implication is that there is direct uptake of atmospheric CO2 (rather than uptake through the water column) by seagrass when seagrass leaves are on the water surface.Atmospheric CO2 is therefore directly taken up within a thin film of water over the seagrass leaves, but this seagrass-driven CO2 flux is not included in the BF flux calculations.
The fact that the CO2 flux was larger than the BF flux may have been partly caused by the temperature and CO2 gradients in the atmospheric layer between the EC measurement height and the water surface.Vertical gradients in air temperature frequently occur because of the large difference in temperature between the atmosphere and water.Indeed, the difference between the temperature at the EC measurement height and the water temperature (ΔT) ranged from +8 °C to −10 °C during the measurement period (Fig. 10).The atmospheric CO2 concentration was inversely Page 18 of 40 Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-499Manuscript under review for journal Biogeosciences Discussion started: 14 December 2017 c Author(s) 2017.CC BY 4.0 License.
proportional to air temperature, mainly because of the ideal gas law.The inflow of air masses with different temperatures and CO2 concentrations from other regions surrounding the lagoon (e.g., terrestrial environments and open sea) could also have contributed to this inverse relationship.We therefore believe that the ΔT caused a vertical gradient in atmospheric CO2 concentration, the result being a vertical CO2 flux in accord with Fick's law.Indeed, the absolute values of the EC CO2 fluxes were large when ΔT was large (Fig. 11b).For example, a large negative EC flux that was more than 10 times the maximum negative BF flux was observed when the ΔT was negative and large in magnitude.This phenomenon frequently occurs when a cold air mass overlies a warm air mass (Fig. 11a).The opposite phenomenon is observed when a warm air mass overlies a cold air mass (Fig. 11b).The positive temperature gradient did not result in a large EC flux similar to the flux associated with the negative temperature gradient because the positive temperature gradient produced stable stratification and thereby prevented vertical eddy movement and CO2 flux, a phenomenon observed in the case of terrestrial EC fluxes at night (Aubinet and Feigenwinter, 2010).Nevertheless, a temperature gradient-driven vertical CO2 flux may explain the discrepancy between the EC and BF fluxes.
In summary, we attribute the discrepancy between the EC and BF fluxes to the following four factors: (1) major technical uncertainties in both methods; (2) differences in measurement location and in the temporal scales of the measurements; (3) limitations of the BF method related to the presence of vegetation on the water surface; and (4) the gradient of atmospheric conditions between the height of the EC measurements and just above the water surface (Table 1).The latter two factors may cause the EC CO2 flux to   Detrending was performed by using block averaging (BA).Our new PP method (PP2) included three data-filtering steps based on the received signal strength indication of the CO2 sensor and the standard deviation of the CO2 and water vapor concentrations divided by the corresponding average standard deviation during the measurement period (nSD).Detrending in PP2 was performed by using high-pass filtering (Massman, 2000).
Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-499Manuscript under review for journal Biogeosciences Discussion started: 14 December 2017 c Author(s) 2017.CC BY 4.0 License.The PP2 procedure combines a series of filtering methods based on the received signal strength indicator (RSSI) of the EC sensor, the normalized standard deviations (nSDs) of the atmospheric CO2 and water vapor concentrations, and high-pass (HP) filtering.

(
http://www.licor.com/env/newsline/2012/04/overcoming-the-challenges-of-open-pathmethane-measurements-with-the-li-7700/).Second, the criteria for excluding erroneous fluxes were identified.Erroneous fluxes were identified based on unnatural spikes, jumps, and shifts in the data.Outliers were excluded based on three statistical parameters: (1) the normalized standard deviation, nSD (i.e. the SD over a 30-min period divided by the average value for the entire measurement period after RSSI filtering; (2) skewness; and (3) the absolute value of kurtosis.Each of these threshold values was determined after the top four outliers for each parameter were excluded; the raw data (CO2 concentration, wind speed) associated with the top four outliers were assumed to be erroneous changes that would not happen naturally.Page 10 of 40 Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-499Manuscript under review for journal Biogeosciences Discussion started: 14 December 2017 c Author(s) 2017.CC BY 4.0 License.
(a)  shows the retained CO2 flux data.The largest positive CO2 flux (release to atmosphere) was 156.51 μmol m −2 s −1 at 2:00 on 23 June (day 56).The largest negative CO2 flux (uptake of atmospheric CO2) was −217.93 μmol m −2 s −1 at 22:00 on 4 October (day 129).These fluxes were more than three orders of magnitude larger than the magnitude of the average of the measured EC fluxes.Figure4shows the instantaneous atmospheric CO2 concentration, water vapor concentration, and the cumulative covariance between CO2 and vertical wind speed during the times when the CO2 fluxes were most positive or most negative.Page 13 of 40 Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-499Manuscript under review for journal Biogeosciences Discussion started: 14 December 2017 c Author(s) 2017.CC BY 4.0 License.
3(b).Of the 2493 total data points remaining after PP1, approximately 234 (9 %) were excluded by RSSI filtering.Subsequent nSD filtering removed 426 additional data points (17 %); approximately 73 % of the measurement data remained after this filtering.The mean and SD of EC CO2 flux after PP2 were −0.54 and 2.2 μmol m −2 s −1 , respectively.The nSD threshold was calculate using the averages of the CO2 and water vapor concentrations for the entire measurement period after RSSI filtering (CO2: 16.02 mmol m −3 , vapor: 548.10 mmol m −3 ).The nSD threshold value obtained from the average values was 0.050.Other thresholds were 0.48 for skewness and 3.1 for the absolute value of kurtosis.Among these comparisons, the nSD was determined to provide the best threshold for the measurement data because the number of data remaining after filtering was the largest for the nSD (nSD: 1661 [74 %], skewness: 422 [19 %], and the absolute value of kurtosis: 445 [20 %]).The nSD filtering resulted in the exclusion of 137 (6 %) of the CO2 data points and 569 (25 %) of the water vapor data points.Figure 5 shows the three parameters for the EC CO2 flux data.HP filtering decreased the absolute values of some CO2 fluxes (Fig. 6).Most of the data were not changed very much by HP filtering; however, HP filtering decreased the Page 14 of 40 Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-499Manuscript under review for journal Biogeosciences Discussion started: 14 December 2017 c Author(s) 2017.CC BY 4.0 License.

Page 16 of 40
Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-499Manuscript under review for journal Biogeosciences Discussion started: 14 December 2017 c Author(s) 2017.CC BY 4.0 License.filterout the noise associated with the spatial heterogeneity of the BF fluxes.In July, the BF fluxes estimated with the W92 seemed to be the most consistent with the EC fluxes.
be larger than the BF flux in aquatic systems that have large amounts of vegetation or Page 19 of 40 Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-499Manuscript under review for journal Biogeosciences Discussion started: 14 December 2017 c Author(s) 2017.CC BY 4.0 License.are located where temperature gradients form in the atmospheric boundary layer.Determination of the contribution of aquatic ecosystems to mitigating the adverse effects of climate change will require consideration of all processes related to atmosphere-aquatic ecosystem exchange.For this purpose, the EC CO2 flux should be a more robust indicator than the BF flux, which includes only processes related to airwater exchanges.Improving the EC method is therefore essential for a re-evaluation of atmosphere-aquatic ecosystem CO2 gas exchanges and comprehensive analyses of the contributions of aquatic environments to mitigating the adverse effects of climate change.Competing interests.The authors declare that they have no conflict of interests.Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-499Manuscript under review for journal Biogeosciences Discussion started: 14 December 2017 c Author(s) 2017.CC BY 4.0 License.Table 1.Summary of the differences in fluxes calculated by the EC and BF methods , https://doi.org/10.5194/bg-2017-499Manuscript under review for journal Biogeosciences Discussion started: 14 December 2017 c Author(s) 2017.CC BY 4.0 License.

Figure captions Figure 1 .
Figure captions

Figure 2 .
Figure 2. Post-processing (PP) methods for EC flux calculation.The EC CO2 flux was calculated by using the conventional PP method (PP1) in EddyPro.The corrections involved in PP1 have been described in previous publications (e.g., Lee et al., 2004).

Figure 4 .
Figure 4. Instantaneous atmospheric CO2 concentration (a), water vapor (atmospheric H2O) concentration (b), and cumulative covariance of atmospheric CO2 concentration and vertical wind speed calculated with PP1 when the CO2 fluxes (c) showed the largest

Figure 5 .
Figure 5.Comparison between the effects of the three filtering parameters and the CO2 flux after the RSSI filtering procedure.The red circles indicate the top four outliers for each parameter, which were determined as erroneous (not natural) fluxes.Each threshold (broken red line) was determined so as to remove these data.(a) Normalized standard deviation (nSD = standard deviation over 30-min divided by the average during the entire measurement period; threshold = 0.05; 74 % of data retained).(b) Skewness (threshold = 0.48; 19 % of data retained).(c) Absolute value of kurtosis (threshold = 3.1; 20 % of data retained).

Figure 6 .
Figure 6.Comparison of CO2 fluxes calculated by the block averaging (BA) method and EC flux with HP filtering.Most of the data lay on or close to the solid line (y = x)and were not much changed by HP filtering.However, in the case of the data in the shaded area, HP filtering decreased the absolute value of fluxes by removing the longterm effect of CO2 change (see Fig.7).

Figure 7 .
Figure 7. Examples of the deviation calculations of atmospheric CO2 concentration (a) and the cumulative covariance of atmospheric CO2 concentration and vertical wind speed (b).

Figure 8 .
Figure 8. Power spectra of CO2 flux with PP1 and PP2.The spectra were normalized

Figure 9 .
Figure 9.Comparison of BF flux, EC flux with PP1, and EC flux with PP2 in May (a), July (b) and September (c).

Figure 10 .
Figure 10.(a) Relationship between air temperatures and atmospheric CO2 concentrations.The gray circles indicate the data averaged over 30 min.The open diamonds and error bars indicate the binned averages every 1 °C and 1 SD, respectively.The solid line indicates the slope estimated from the change in air volume assuming that CO2 behaves as an ideal gas.(b) Relationship and between temperature difference (ΔT) and CO2 flux.The gray circles indicate the EC flux data calculated every 30 min.The solid diamonds and error bars indicate the binned averages every 1 °C and 1 SD, respectively.Note that for clarity not all data are plotted on both graphs.

Figure 11 .
Figure 11.Schematic diagram showing the relationships between negative temperature gradient and negative EC CO2 flux (a) and between positive temperature gradient and positive EC CO2 flux (b).The atmospheric CO2 gradient reflects the air temperature gradient and assumes that CO2 behaves like an ideal gas.The temperature gradientdriven vertical CO2 flux may explain the discrepancy between the EC and BF fluxes.