This paper reports a fairly long and potentially valuable time series of CO₂ and CH₄ fluxes measured directly by the eddy covariance (EC) method over a boreal river. EC method certainly offers many advantages over other methods (chamber, tracer) in terms of temporal and spatial scales. At the same time, careful processing and filtering of EC data are critical for the data to be trust worthy. I have some reservations about the authors’ flux processing/filter methods. Please see detailed comments below.

Concurrent observations of pCO_2w, water flow speed and discharge, in addition to winds, seem to offer the ideal setup for examining the drivers for the transfer velocity of CO₂ (flux/dC). Yet puzzlingly, the authors have chosen to bypass this calculation entirely. The multi-linear regression approach has some usefulness for investigating the drivers of pCO_2w. But it’s less suitable for investigating the drivers of flux in my opinion. As it stands, I don't feel like I've learned much new from this paper.

That the CO2 flux measurements were made from an undried Licor raises question marks, as it’s well known by now that H₂O can cause a substantial interference in CO₂ flux over water. If the Picarro also measured CO₂, the authors can compare CO₂ flux from the two instruments.

I have a strong suspicion that the flux footprint overlapped onto the land upwind at times, which possibly contributed to the observed diurnal cycle in CO₂ flux. Please check carefully.

• We thank Dr. Yang for his insightful and detailed review. We address these questions as well as all the detailed comments below. Please see our response below.

Detailed comments:

Line 44-46. Not sure if this brief section warrants mentioning at this stage.

• We have reworded the latter part of the sentence to make it less bold of a claim.

Line 106. Lack of drying is a concern for the Licor7200. Several works have now shown that for air-sea flux measurements, not drying the air will result in biased CO₂ flux measurements due to imperfect H₂O correction within the Licor (e.g. www.atmos-chem-phys.net/14/3361/2014/ ).

The H₂O issue may be less severe for the Picarro (e.g. www.atmos-meas-tech.net/9/5509/2016/), perhaps due to the very sharp laser frequency (compared to Licor’s broadband).

• According to Nilsson et al. (2018) (https://doi.org/10.1175/JTECH-D-17-0072.1) who studied the effect of H₂O on CO₂ fluxes measured with Li-7200, the effect was largely due to sea salt contamination on the optical windows of the gas analyser. A similar setup without a dryer has been used in other freshwater studies, e.g. by Huotari et al. (2013), Erkkilä et al. (2018) and Spank et al. (2020) (references in the manuscript). It is a standard setup that follows the ICOS protocol (Sabbatini et al. 2018) (https://doi.org/10.1515/intag-2017-0043) and the data processing is well-established. The volume effect of H₂O on CO₂ and CH₄ fluxes has been taken care of by utilising dry mole fractions of CO₂ and CH₄.

Does the Picarro not measure CO₂ as well?  It would be insightful to compare CO₂ flux from the Picarro vs. the Licor7200.

• Our Picarro did not measure CO₂. It would indeed have revealed important insight of the ecosystem and the measurement setup.

Line 122: lag of 0.34 s was used for CO₂ and 7.0 s for CH₄

Looking at the indicated inlet dimensions and flow rates, I compute a delay time of <0.1 s for CO₂ and <0.6 s for CH₄.  These look to be very different from the estimated delay time from the maximum covariance method, especially for CH₄.  Did the Picarro subsample from the 10 m inlet?  Or were too many outlier delay estimates included in the mean calculation?  The computed flux will be biased low if the delay time is incorrect.

• The lag time calculation with the maximum cross-covariance method was limited to a certain time window and there were no significant outliers in the lag times. We recalculated the lag times, and the results were the same as in the manuscript. The lag times of both gas analysers are increased a little from the theoretical value by the sample line filters.  One reason between the discrepancy in the case of Picarro could be clock synchronisation between the sonic anemometer and the Picarro. We added a figure in the manuscript showing examples of the cross-covariance functions for both CO₂ and CH₄. The figure is also below.
• 
• Figure 1. Example of the time lag determination for (a) CO₂ and (b) CH₄. The solid line shows the cross-covariance function of the vertical wind w and the gas mixing ratio χ. The vertical dashed lines mark the windows within which the lag times were searched. In this example, the location of the cross-covariance peak, i.e. the time lag for these averaging periods, was 0.3 s for CO₂ and 7.1 s for CH₄.

Line 132. I don’t see this low frequency ‘attenuation’ correction being applied very often. I understand that by linear detrending, an implicit decision is made such that variability below a certain frequency is not considered to be flux in the EC calculation. A typical averaging time scale of 30 min is never going to statistically capture those very slow varying eddies. But this low frequency signal (e.g. due to convection within the boundary layer) could have either a positive or negative contribution of the flux, right? Then in that case is it really an ‘attenuation’?

• The low-frequency correction is a standard procedure and implemented according to the ICOS protocol (Sabbatini et al. 2018). In the form that it is implemented, it aims to correct the flux attenuation from excluding the largest eddies. However, for our site the correction is very small (1–2%) because of the rather low measurement height and the absence of strong convection (on average z<<-L, where L is the Obukhov length).

More generally, it’d be useful to see the mean CO₂ and CH₄ cospectra.

• They are now added to the manuscript.

Line 158. ‘removed’ or ‘retained’?  I would’ve thought that the low variability data should be retained.

• The correct word is ‘retained’. The error is now fixed.

Line 161. z₀ of 0.01 m seems large for a water body. Furthermore, the flux footprint is very sensitive to the stability of the atmosphere. Stable conditions lead to much larger footprints and in this case more overlap with upwind land. From the measured sonic flux, the authors can probably screen out cases when the atmosphere is stable. I’m guessing that for a lot of the daytime measurements, T_a > T_surf and it was fairly stable.

• We had another eddy covariance system on the platform measuring fluxes at 1 m height adjacent to the system described in the manuscript. Calculating z₀ from the two wind measurements resulted in z₀ ≈ 0.01 m. The straight section in the river and the fetch being close to or exceeding 1 km enabled rather large wind-generated waves, especially when the wind was from the south, i.e. against the current. At highest, the waves were approximately 0.5 m tall and 0.1–0.2 m waves were common. In this sense, z₀ = 0.01 m is not unreasonable.
• In footprint calculations, all meteorological characteristics are considered. The acceptable wind directions regarding the fetch are following the calculation results. In addition, we desired to be conservative and only directions which are more along the river were selected, although footprints would have indicated acceptable fetch for directions more towards the bank. Note also that the footprint probability density function is weighted near the measurement location. Thus, even stable cases allow for acceptable data as the measurement height was small.

Fig 6a. that CO₂ flux goes negative in midday in June, to me, suggests that the EC flux footprint is overlapping onto land further upwind. (As shown in Figure 3, pCO_2w >> pCO_2a during the entire campaign). This is probably because of the stable atmosphere when T_a > T_surf. I suggest the authors check their flux filtering more carefully before commenting further on the apparent diurnal cycle in CO₂ flux. 

• See the previous answer. Even though we intended to be careful in data filtering regarding the fetch, it is true that some flux signal – beyond the footprint extent of 90% – originates from land. This creates scatter in the data but it cannot change the sign of the flux. The mismatch between the flux sign and pCO_2w/pCO₂ ratio may result from several reasons, like that the concentration measurement is not fully representing the larger flux footprint.

Section 4.1 see this paper on filtering of CO₂ and CH₄ fluxes from a coastal site, which may be of use: www.atmos-chem-phys.net/16/5745/2016/

In particular, the authors could use dC/dt and <u’C’> as additional filtering criteria for longer distance transport.

• We have closely followed the ICOS protocol (Sabbatini et al., 2018) in filtering the fluxes. Additionally, we have utilised the approach of further calculating the mean wind direction in 5-minute intervals and using those wind directions as additional filtering criteria (Section 2.3, lines 166–179). According to Blomquist et al. (2012) (https://doi.org/10.5194/amt-5-3069-2012) who used dC/dt and <u’C’> as criteria for filtering fluxes and who are cited in the paper above, these criteria are used for ensuring the stationarity in the time series. However, we have already utilised the flux stationarity criterion by Foken and Wichura (1996). dC/dt and <u’C’> might be suitable for assessing the effect of lateral advection in the fluxes, but we have implemented the gas mixing ratio standard deviation filtering for excluding cases with strong advection.

Also, I’m not sure if section 4.1 is most suitable for discussion. Might be more suitable for the methods section.

• This section contains discussion about applying the eddy covariance method over water bodies rather than describing the method itself. For example, the section contains discussion about the 5-minute wind direction method. To expand the discussion, we have added a reference to the paper by Blomquist et al. (2012) to Section 4.1.

Table 5. I understand that pCO_2a didn’t vary hugely and thus variability in dpCO₂ is more driven by variability in pCO_2w. Still, shouldn’t be the air-water concentration difference in pCO₂, rather than pCO_2w that’s most relevant here?

• This is a valid point. The multivariate model calculations were redone using dpCO₂ instead of pCO_2w. The resulting models and their ranking are almost identical to the ones in the preprint. The R² and p values changed only slightly. The tables in the manuscript are now updated accordingly. The emergent picture of the flux drivers did not change and thus the discussion remains the same.

Also, why don’t the authors compute the CO₂ transfer velocity (flux/d_C), and then explore how the transfer velocity may be affected by these secondary terms? It seems strange to have multi-linear relationships where flux is parametrized as e.g.: a * pCO_2w + b*U.  Surely something like flux = (a + b * U + c*U_w) * dpCO_2w would be better (I’ve neglected temperature dependencies for simplicity).

• We are writing another manuscript related to the gas transfer velocity using the same dataset and where we will analyse in detail the dependence of the transfer velocity on e.g. flow speed. The purpose of our multivariate analysis in this manuscript with daily and weekly means is to gain a basic understanding of flux drivers.

Line 351. Is there a diurnal cycle in measured pCO_2w? That’d help to address this inconsistency.

• There is a small diurnal cycle in the measured pCO_2w. The maximum mean difference in pCO_2w is 46 ppm. A figure is of the mean diurnal pCO_2w is now added in the manuscript.

Line 364. Were there previous pCH_4w measurements at this site?  Or pCH₄ from similar rivers?  What sort of CH₄ supersaturation and flux are we expecting roughly?

• There are no previous measurements of pCH_4w at this site. We have added a paragraph in section 4.2 where we summarise earlier studies of pCH_4w and FCH₄ from similar rivers.

The preprint provides an interesting dataset of eddy covariance flux measurements at a river site in Finland and as the paper states eddy covariance measurements of CO₂ and CH₄ at rivers have only been done once before. Such measurements are especially difficult due to the inhomogeneous surroundings at small to medium rivers and have to be analysed carefully.

• We thank Dr. Zavarsky for his thoughtful and thorough remarks. Please see our responses below.

Scientific questions:

Line 40: “The River Kitinen is a regulated river and, as such, provides information on e.g. how the anthropogenic modification of the flow affects gas exchange”

 

Have you extensively addressed this later in the paper? Later there is only a short sentence about this.

• We decided to remove this sentence from this manuscript and will address this question in the second manuscript focusing on the gas transfer coefficient.

Line 159: You are using the Kljun et al 2015 footprint calculations. Did you put in the whole time series of wind speed, stability etc and the calculate the footprint climatology? Did you exclude certain “good directions” (parallel to the river) EC data because of an extensive footprint? What does “the river bank directions are excluded” mean? Did you exclude them anyways or was the footprint calculation stopped there?

• The full measured dataset of the mentioned variables was used in the calculation and periods when the wind was from the river bank sectors are excluded from further flux analysis. We rephrased this paragraph to better explain the calculation process. We want to add here that the data filtering is a compromise between data statistics and the safest fetch. Narrowing even more the acceptable wind directions would result into a very low amount of data points. See the comments above related to Dr. Yang’s remarks about the footprint.

Line 253: You chose two timespans for day and night-time. Do you have a total radiation data set to back this up or is it a feature of your EC data set? You measured PAR in the water but some questions (also by Mr Yang) refer to land vegetation. A total radiation dataset would be helpful. Wasn’t there a meteorological mast nearby? Was there such a sensor mounted?

• This rudimentary classification was selected to get an idea of diurnal differences. We are not using the classification for any further analysis. We indeed have a full radiation dataset, but for avoiding land contamination in the fluxes, we have already implemented several methods: flux stationarity criterion (Foken and Wichura, 1996), limiting the standard deviations of the gas mixing ratios and the new method of observing the 5-minute mean wind directions.

You have correlated fluxes, PAR and, pCO_2w on a daily scale. What about hourly/2 or 3 hourly scale fluxes with respect to PAR or total radiation? This could also help interpreting your diel cycles.

• We have conducted an additional multivariate analysis on average diel cycles of the same variables as in the daily and weekly analyses in the manuscript (excluding precipitation, as we only had daily precipitation rates) and will add those results and discussion in the manuscript. According to this analysis, pCO₂ correlates the most with T_surf and U_w. For FCO₂, the main drivers are pCO₂ and U and for FCH₄, they are T_surf and pCO₂.

Line 327: Good point. Sections across the river with pCO_2w values would be interesting.

• We agree, however, we were measuring pCO_2w in only one location during the campaign.

Line 383: there a plenty of gridded precipitation reanalysis products. One could combine that with watershed geographical data and get a better picture of the upstream precipitation.

• This would require a modelling approach and extensive research on the river discharges in the watershed, flushing of carbon from soils etc. and is beyond the scope of this manuscript. Our approach focuses instead on simple environmental variables that are relatively easy to obtain.

Technical Questions:

Line 227: “The 30-minute average air temperature rose from 10-20°C …. to 30°C. “Didn’t it just rise from 10°C to 30°C.

• This sentence is reworded to make it clearer to the reader.

Figure 3:  I would recommend changing the water temperature heat plot to a line plot with maybe surface, depth at x m and bottom. I couldn’t figure out any temperature features over depth. It is mostly variable over time.

• We have changed the heat plot to a line plot.

Additionally, a total radiation dataset, as mentioned above, would be helpful.

• We would like to show in the manuscript only those variables that are used in the analysis but a full radiation dataset can be added in the supplementary information.

Line 269: In this paragraph you state the you present the five best models. You could make them more visible and distinguishable: like “a)….b)…c)…”

• The letters are now added to Tables 4, 5 and 6.

Line 273: “Some of best daily models include U”. Can you be more precise.

• We now mention the exact number of models of the five best ones than contain U as a driving variable.

Line 366: This paragraph covers a lot. Maybe some subsections or partitions would be helpful.

• We now divide this section into drivers of pCO₂ and drivers of fluxes to make it more readable.

Line 397: “…on only some of the best models”. Can you be more precise.

• Similarly as in the comment above, we now mention the exact number of models containing U.