Measuring turbulent CO₂ fluxes with a closed-path gas analyzer in a marine environment

In this study, we introduce new observations of sea–air fluxes of carbon dioxide using the eddy covariance method. The measurements took place at the Utö Atmospheric and Marine Research Station on the island of Utö in the Baltic Sea in July–October 2017. The flux measurement system is based on a closed-path infrared gas analyzer (LI-7000, LI-COR) requiring only occasional maintenance, making the station capable of continuous monitoring. However, such infrared gas analyzers are prone to significant water vapor interference in a marine environment, where CO₂ fluxes are small.

Two LI-7000 analyzers were run in parallel to test the effect of a sample air drier which dampens water vapor fluctuations and a virtual impactor, included to remove liquid sea spray, both of which were attached to the sample air tubing of one of the analyzers. The systems showed closely similar (R²=0.99) sea–air CO₂ fluxes when the latent heat flux was low, which proved that neither the drier nor the virtual impactor perturbed the CO₂ flux measurement. However, the undried measurement had a positive bias that increased with increasing latent heat flux, suggesting water vapor interference.

For both systems, cospectral densities between vertical wind speed and CO₂ molar fraction were distributed within the expected frequency range, with a moderate attenuation of high-frequency fluctuations. While the setup equipped with a drier and a virtual impactor generated a slightly higher flux loss, we opt for this alternative for its reduced water vapor cross-sensitivity and better protection against sea spray. The integral turbulence characteristics were found to agree with the universal stability dependence observed over land. Nonstationary conditions caused unphysical results, which resulted in a high percentage (65 %) of discarded measurements. After removing the nonstationary cases, the direction of the sea–air CO₂ fluxes was in good accordance with independently measured CO₂ partial pressure difference between the sea and the atmosphere. Atmospheric CO₂ concentration changes larger than 2 ppm during a 30 min averaging period were found to be associated with the nonstationarity of CO₂ fluxes.

1 Introduction

Anthropogenic actions, such as combustion of fossil fuels and land use changes, have perturbed the global carbon cycle, resulting in climatic changes (Solomon et al., 2009). A quarter of the anthropogenic carbon dioxide (CO₂) emissions to the atmosphere are bound by the oceans (Heinze et al., 2015), causing ocean acidification (Feely et al., 2009). To better understand the global carbon cycle, measurements of the CO₂ exchange between the atmosphere and marine ecosystems are required. These measurements are also useful for developing parameterizations of gas exchange intensity, such as the gas transfer velocity, used in global carbon models (Takahashi et al., 2002). In the case of coastal seas, sea–air CO₂ flux measurements provide information about the feedbacks between the aquatic carbonate system and marine ecosystems, since the direction and magnitude of these fluxes depend on the photosynthetic carbon assimilation in surface seawater. Although the continental margin seas cover only a small portion of the oceans, up to 15 % of the total ocean primary production takes place in these seas, which are responsible for over 40 % of the total oceanic carbon sequestration (Muller-Karger et al., 2005).

With the development of fast-response infrared CO₂ analyzers, it has become possible to apply the eddy covariance method for directly measuring the sea–air CO₂ fluxes (Jones and Smith, 1977). Early on, however, Webb et al. (1980) showed that the CO₂ flux measurements made with this technique need to be corrected for the effects of temperature and water vapor (H₂O) fluctuations. This correction can be significant in aquatic environment, even on the same order as the measured CO₂ flux (Sahlée et al., 2011). More recently, it has been recognized that CO₂ gas analyzers suffer from water vapor cross-sensitivity (Kohsiek, 2000). Blomquist et al. (2014) concluded that this is the most significant error source in sea–air CO₂ flux measurements, because the CO₂ fluxes in a marine environment are typically small, as compared to terrestrial ecosystems. A solution to the cross-sensitivity problem is to dry the sample air before the measurement (Miller et al., 2010). Also, efforts have been made to correct for the cross-sensitivity problem in the data post-processing phase (Edson et al., 2011; Prytherch et al., 2010). However, these corrections may be inadequate and thus do not obviate the use of a drier (Landwehr et al., 2014).

The cross-sensitivity results from the overlap of the infrared radiation frequency bands of CO₂ and H₂O, because of which the H₂O molecules present will increase the apparent CO₂ concentration. If the optical filter, used for the selection of the transmitted frequency band, leaks out-of-band radiation, a substance with a different absorption frequency band can interfere with the measurement. This effect is referred to as the direct absorption cross-sensitivity. In addition, the collisions with different molecules cause the frequency bands to broaden (so-called pressure broadening). By testing some commercial infrared gas analyzers, Kondo et al. (2014) found that the factory-calibrated correction for the direct absorption interference may not be optimized and that the pressure-broadening effect caused an overestimation of the CO₂ flux, which increased with increasing water vapor flux.

Infrared gas analyzers are classified according to the type of the optical path: the open-path analyzers measure the absorption of the infrared signal in ambient air, whereas the closed-path gas analyzers have an enclosed measurement chamber. Both types suffer from H₂O cross-sensitivity but otherwise have differing pros and cons; for instance, the closed-path instruments are known to act as low-pass filters, which generate a loss in the measured flux (Leuning and King, 1992). On the other hand, a long sample line attenuates temperature fluctuations, thus eliminating the need for correcting for sample air expansion/contraction (Rannik et al., 1997). For both analyzer types, the dilution due to water vapor (Webb et al., 1980) can be corrected accurately as a point-by-point operation on the high-frequency data (Ibrom et al., 2007). However, such an approach is not feasible for compensating the temperature fluctuations with open-path sensors because of the spatial separation between the CO₂ and the sonic temperature measurements.

As direct sea–air gas exchange measurements are performed in the surface boundary layer, in a close proximity to the water surface, liquid sea spray may block the optical path of open-path sensors and clog the inlet of closed-path analyzers. In addition, high relative humidity (RH) can produce H₂O condensation on the lenses of an open-path system. In a closed-path system, the accumulation of sea salt on the optical lenses is minimized as the inlet is typically protected by a Teflon filter. Efforts have been made to solve the sea spray contamination problem of open-path sensors by cleaning the optics regularly (Kondo and Tsukamoto, 2007). However, this may be technically challenging, because sea salt films can be formed on the windows of an open-path analyzer in a matter of hours (Miller et al., 2010). Sample air drying has been recommended as a solution to the sea salt contamination (Nilsson et al., 2018).

Open-path gas analyzers have been mostly applied for sea–air CO₂ flux measurements due to their low power consumption, small high-frequency attenuation and ease of data synchronization with wind measurements (Blomquist et al., 2014). Additionally, Miller et al. (2010) found that closed-path sensors are more sensitive to (ship) motion than the open-path sensors. However, this does not apply to all closed-path analyzers, since some of them have internal design resembling open-path instruments (LI-COR, 2016).

While the eddy covariance method is widely used for directly measuring the surface–atmosphere exchange of energy and matter, it is based on several theoretical assumptions. These include the horizontal homogeneity of terrain and the stationarity of transport processes, and that turbulence is fully developed and there exist no transport mechanisms other than vertical turbulence (Dabberdt et al., 1993). Foken and Wichura (1996) noted that flow nonstationarity is one of the most serious problems affecting the surface exchange measurements, as in such conditions the observed turbulent flux does not equal the flux at the surface. The stationarity assumption can be violated owing to diurnal forcings and varying weather patterns, for example. Horizontal CO₂ gradients develop easily in the vicinity of the sea–land boundary, violating the assumptions of horizontal homogeneity and negligible advection (Rey-Sánchez et al., 2017). Because an obstacle in a flow field generates an uplift zone on the windward side of an obstacle (Chatziefstratiou et al., 2014), there should not be steep elevation changes near the flux tower in order to avoid nonzero vertical wind speed. Also the roughness and flux discontinuities on a sea–land boundary generate secondary circulation that drives advection and causes challenges for coordinate rotation (Kenny et al., 2017).

The Baltic Sea forms a large and diverse biogeochemical system, providing a great potential for studying interactions between the marine ecosystem and the aquatic carbonate system. However, only a few fixed measurement sites measuring sea–air CO₂ fluxes are located in the Baltic Sea (Lammert and Ament, 2015; Rutgersson et al., 2008). A micrometeorological tower placed on the shore of an island offers a cost-effective alternative to ship measurements, as maintenance and installations are more feasible. Moreover, a motion correction, required with ship measurements, is not needed, and the flow distortion can be minimized with a suitable positioning of the flux tower.

In this paper, we introduce a newly established eddy covariance measurement site, located on the island of Utö in the Baltic Sea. This paper has two objectives: (1) to study the characteristics of the new site and measurement setup and (2) to analyze empirically the effect of water vapor on the CO₂ sea–air fluxes. After experimenting with an open-path and a semi-open-path gas analyzer, we opted for a closed-path sensor, which can be protected against sea spray contamination and allows straightforward implementation of sample air drying. We compare measurements made with two identical closed-path analyzers, one of which is equipped with a drier and a virtual impactor. We also address the quality of these measurements by analyzing the stationarity and integral turbulence characteristics of the flow and the homogeneity of flux footprints.

2 Materials and methods

2.1 Site description

The measurement site is located on the island of Utö on the southern edge of the Archipelago Sea in the Baltic Sea (59^∘46^′55${}^{\prime \prime }$ N, 21^∘21^′27${}^{\prime \prime }$ E) (Fig. 1). The Archipelago Sea is a small sea area between the southwest coast of Finland and the sea of Åland and is characterized by thousands of small islands. The island of Utö is a treeless cliff with small shrubs, with an area of 0.81 km². Since 2012, the Finnish Meteorological Institute has been building up the new Utö Atmospheric and Marine Research Station in collaboration with the Finnish Environment Institute (Laakso et al., 2018). The choice of this island was due to its easy access, well-developed technical infrastructure and permanent inhabitation, as well as the existing long-term meteorological, marine and air quality measurements.

Figure 1 Location of the island of Utö in the Archipelago Sea. The research installations on the island consist of the Utö Atmospheric and Marine Research Station and its flux tower (A), the inlet of the seawater sampling system (B), the Atmospheric ICOS station (C), and a weather and air quality station (D). The red solid lines indicate the wind sector used for sea–air flux measurements (180–260^∘) used in this paper, and the dashed line indicates the fetch-limited sea sector excluded from the current analysis.

The Marine Research Station and its flux tower are located on the western side of the island. To the west of the shore, the water depth quickly deepens to 80 m. The closest shoal, Tratten, is located 1.4 km to the southwest and has an area of 0.007 km².

The mean annual wind speed at Utö is 7.1 m s⁻¹ (Pirinen et al., 2012). On average, the minimum monthly wind speed occurs in July (5.6 m s⁻¹) and the maximum in December (8.9 m s⁻¹). The wind blows approximately 60 % of the time from the sector between south and northwest (Fig. 1), providing a good temporal coverage for the sea–air flux measurements.

The mean annual air temperature in Utö is 6.6 ^∘C, whereas the mean annual surface seawater temperature at the depth of 5 m is 8.1 ^∘C (Laakso et al., 2018). In deep water (at 50 m depth), the mean annual seawater temperature is 4.4 ^∘C. The maximum monthly mean air temperature (16.7 ^∘C) is observed in July, whereas the seawater at 5 m depth reaches its maximum temperature (16.1 ^∘C) in August. The minimum temperature of air (−2.1 ^∘C) occurs in February, while the temperature of the surface seawater is at its minimum (1.0 ^∘C) in March. On average, there is ice around Utö every few years, with a typical ice cover duration ranging from 1 to 3 months.

Unlike the oceans, the carbon cycle in the Baltic Sea is heavily influenced by biological activity. The partial pressure of CO₂ (pCO₂) in surface water has its maximum 60 Pa (600 µatm) during winter, when biological activity is diminished, mineralization prevails and mixing processes transport CO₂-rich seawater to the surface (Wesslander, 2011). During summer, pCO₂ in surface water declines to approximately 10 Pa (100 µatm) as a result of primary production. Due to the annual cycle of primary production, the Baltic Sea is a source of CO₂ to the atmosphere in winter and a sink in summer. Upwelling and the diurnal biological cycle generate short-term variations in the surface water pCO₂.

2.2 Instrumentation

A 9 m tall micrometeorological tower was placed on the western shore of the island (Figs. 1 and 2). The tower is mounted on a cliff that has no vegetation. The base of the tower is approximately 3 m above the mean sea level, and the horizontal distance between the tower and the sea is approximately 4 m. In this paper, we utilize data from the open sea sector towards the Baltic Proper (180–260^∘). The fetch in the sector 260–340^∘ may be limited, and, as the anemometer base potentially disturbs the flow, the data from this sector were not included in this study.

Figure 2 The flux tower (height 9 m) on the shore of the island of Utö during a wintertime westerly storm (photo by Ismo Willström).

Air velocity components together with air temperature, T, were measured with an acoustic anemometer/thermometer (uSonic-3, Metek) attached to the top of the tower (11.5 m a.s.l.). CO₂ and H₂O molar fractions were measured by using two closed-path differential non-dispersive infrared gas analyzers (LI-7000, LI-COR) that were placed in a 1.5 m tall instrument hut close to the tower (Fig. 2). The inlets of both sample flows are protected with a grate and are located directly beneath the anemometer; the distance between the inlets and the lower anemometer transducers is 30 cm. The outside parts of the tubes are approximately 10 m long and are made of Teflon (inner diameter of 3.175 mm) and steel (inner diameter of 4.0 mm). The steel-tube sample line is equipped with a virtual impactor (see Appendix A), to protect the instrument from possible exposure to liquid water, and a 30 cm long PD-100T-12-MKA drier (Perma Pure) to attenuate water vapor fluctuations. The drier is based on the partial pressure difference that drives water vapor from the sample air to the purge stream through a Nafion membrane. The sample air leaving the gas analyzer of the test setup is directed to the purge cell of the drier. Thus, the total pressure difference drives the attenuation of water vapor fluctuations. Nafion driers are found to have only a small permeability to CO₂ (Welp et al., 2013).

Inside the instrument hut, the flow continues in 1 m long Bev-A-Line tubes, which are protected by Teflon filters and connected to the analyzers. Data are logged at a 10 Hz frequency and transmitted through the Internet to a server. The hut temperature is regulated using a fan and a radiator, since large temperature changes influence the performance of the CO₂ analyzers.

An external pump is used for producing a sample flow at a rate of approximately 6 L min⁻¹ for both setups. We operate the LI-7000 gas analyzers in the absolute mode, which means that a small stream of zero gas (0 ppm CO₂) is constantly flowing through the reference cell of both analyzers. The same zero gas flows through both analyzers. The actual zero gas bottle (40 L) is located at the station, and a small stream of zero gas flows continuously to the gas analyzers in the instrument hut. The gas analyzers are calibrated with zero and span (364.4 ppm CO₂) gases every 3 months. After each calibration, the Teflon filters are renewed and the virtual impactor is cleaned. The outside tubes are washed with an isopropyl alcohol–water mixture annually. At the same time, the Bev-A-Line tubes and the internal chemicals of the gas analyzers are renewed, and the optical paths are cleaned.

The setup with the Nafion drier and the virtual impactor is referred here to as the “test setup”, and the other setup as the “standard setup”. The standard setup represents a typical closed-path gas analyzer-based eddy covariance setup commonly applied in terrestrial ecosystems (Aurela et al., 2015). The test setup is an improved version of this configuration, designed especially for a marine environment where water vapor and sea spray are likely to give rise to problems.

2.3 Data processing

2.3.1 Flux calculation

As the LI-7000 gas analyzer does not directly provide molar fractions with respect to dry air, we first calculated the dilution-corrected CO₂ molar fraction as

$$\begin{array}{}\text{(1)}& c={\displaystyle \frac{c_w}{\mathrm{1}-q}},\end{array}$$

where c_w is the uncorrected molar fraction of CO₂ and q is the molar fraction of H₂O.

It is assumed that the long pipe lines, both Teflon and steel, strongly attenuate the temperature fluctuations, and thus there is no need for correcting for expansion/contraction effects in the sample air (Rannik et al., 1997). The eddy covariance mass flux was calculated from the fluctuations of molar fractions of an atmospheric quantity, χ, and vertical wind speed, w:

$$\begin{array}{}\text{(2)}& F={\mathit{\rho }}_{\mathrm{d}}\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\chi }}^{\prime }},\end{array}$$

where ρ_d is the dry-air density. The primes denote the 10 Hz fluctuations with respect to a time average, and the overbar denotes arithmetic averaging. An averaging period of 30 min was adopted. For each period, the time delay between the acoustic anemometer and gas analyzer data flows was determined by maximizing the (absolute) covariance in Eq. (2) within a predefined time window, and a double coordinate rotation was applied to ensure zero vertical wind speed.

The fulfillment of several theoretical requirements for the eddy covariance measurements was considered. The nonstationarity of the flux measurements was analyzed by comparing the average of the CO₂ fluxes calculated for 5 min subperiods ($\frac{\mathrm{1}}{\mathrm{6}}{\sum }_{i=\mathrm{1}}^{\mathrm{6}}{\stackrel{\mathrm{‾}}{w^{\prime }{c}^{\prime }}}_{\mathrm{5}\min ,i}$, where i is the index of the subperiod) with the 30 min flux (${\stackrel{\mathrm{‾}}{w^{\prime }{c}^{\prime }}}_{\mathrm{30}\min }$) similarly to Foken and Wichura (1996). The relative nonstationarity, RN, for a given variable, ζ, was defined as

$$\begin{array}{}\text{(3)}& {\text{RN}}_{\mathit{\zeta }}=\left|{\displaystyle \frac{\frac{\mathrm{1}}{\mathrm{6}}{\sum }_{i=\mathrm{1}}^{\mathrm{6}}{\mathit{\zeta }}_{\mathrm{5}\min ,i}-{\mathit{\zeta }}_{\mathrm{30}\min }}{{\mathit{\zeta }}_{\mathrm{30}\min }}}\right|.\end{array}$$

In the case of CO₂ flux, $\mathit{\zeta }=\stackrel{\mathrm{‾}}{w^{\prime }{c}^{\prime }}$. Here, a 30 min period is regarded as stationary if both of the CO₂ flux measurements fulfill the condition ${\text{RN}}_{\stackrel{\mathrm{‾}}{w^{\prime }{c}^{\prime }}}<\mathrm{0.3}$. To test the development of turbulence, integral turbulent characteristics of temperature and vertical velocity were compared to universal surface-layer functions (Thomas and Foken, 2002). It is assumed that turbulent transport of all scalars (T, CO₂ and H₂O) are similar, and thus only the characteristics of T are examined here. For this examination, only stationary conditions within the sea sector are included. An observation is also discarded if a positive momentum flux is measured, as these can be considered a sign of flow nonstationarity (Yang et al., 2016).

2.3.2 Spectral analysis

Spectral analysis was carried out in order to evaluate the high-frequency attenuation of the fluxes. High-frequency attenuation was examined by calculating the cospectra between w and c (${\text{Co}}_{w{\mathrm{CO}}_{\mathrm{2}}}$) and q (${\text{Co}}_{w{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$) and comparing these with the corresponding cospectra between w and T (Co_wT). For this purpose, it is assumed that the high-frequency end of Co_wT is unattenuated and that the atmospheric turbulent transport of heat, water vapor and CO₂ is similar.

Before calculating the half-hourly cospectra using the discrete Fourier transform, the 10 Hz data were linearly detrended and the Hamming window was applied. The cospectra were normalized with the corresponding covariances and interpolated into 64 logarithmically equally spaced bins.

Observations with an appropriate wind direction (180–260^∘) and stationary conditions (${\text{RN}}_{\stackrel{\mathrm{‾}}{w^{\prime }{c}^{\prime }}}<\mathrm{0.3}$) were accepted for spectral analysis. A total of 607 half-hour periods met these criteria. Next, we discarded small fluxes (< 0.1 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), after which we had 381 half-hour periods. Furthermore, a cospectrum was discarded if its peak was not within a predefined range (0.01–0.5 Hz). A total of 283 cospectra passed this frequency range criterion. Also, a cospectrum was discarded if its normalized peak size (${\text{Co}}_{w{\mathrm{CO}}_{\mathrm{2}}}/\stackrel{\mathrm{‾}}{w^{\prime }{\text{CO}}_{\mathrm{2}}^{\prime }}$) was not within predefined range (0.1–1.0). Thus, 238 observations from a 4-month period were used for the spectral analysis.

The spectral attenuation of the flux as a function of frequency, f, was described by an exponential transfer function, Γ:

$$\begin{array}{}\text{(4)}& \mathrm{\Gamma }\left(f\right)=\exp [-\log (\mathrm{2}\left){\left({\displaystyle \frac{f}{f_{\mathrm{0}}}}\right)}^{\mathrm{2}}\right],\end{array}$$

where f₀ is the half-power frequency at which the ratio of cospectral densities is 0.5. Using a nonlinear least squares fit, the transfer function in Eq. (4) was fitted to the ratio of ${\text{Co}}_{w{\mathrm{CO}}_{\mathrm{2}}}$ or ${\text{Co}}_{w{\mathrm{H}}_{\mathrm{2}}\text{O}}$ to Co_wT. To fit the transfer function, the cospectra were normalized so that their peaks were leveled. Also, four outermost points from the high-frequency end were discarded because these points may easily be biased due to the division by very small number and 20 outermost points from the low-frequency end were discarded because these points do not play role in the high-frequency attenuation. In the case of H₂O fluxes, a power-law relationship between f₀ and RH was additionally derived similarly to Mammarella et al. (2009).

To correct for the high-frequency attenuation of fluxes, the universal Co_wT equations reformulated by Horst (1997) from those originally presented by Kaimal et al. (1972) and Kaimal and Finnigan (1994) were used as a reference. As the shape of this universal cospectrum depends on stability and wind speed, the attenuation for each 30 min flux was calculated as a function of these meteorological parameters.

2.4 Auxiliary data

In addition to the flux measurements, we present here data of the dissolved CO₂ concentration in surface seawater, which is measured using a SuperCO₂ system (Sunburst Sensors LLC) connected to a submersible pump. The SuperCO₂ instrument consists of a showerhead equilibration chamber and an infrared gas analyzer (840A, LI-COR). The measured CO₂ molar fraction is corrected for H₂O dilution effect using the measured H₂O molar fraction. Since this CO₂ measurement is not carried out using pre-dried sample air, a H₂O cross-sensitivity effect on CO₂ is possible but likely small compared to the other error sources, such as spatial heterogeneity in seawater CO₂ concentration. To account for sensor drift, the flows of four reference gases are directed into the instrument every 4 h, and a linear correction is calculated from these calibration measurements. The system automatically cleans the equilibration chamber with hydrogen peroxide periodically. The bottom-moored, floating inlet for seawater is located 250 m west of the station (Fig. 1) and is kept approximately 4.5 m below the mean sea level: the sampling depth typically varies within 4.0–5.0 m. The water depth at this location is 23 m. Since sample water exchanges heat with the pipe and its surroundings, a correction for this temperature effect was applied according to Takahashi et al. (1993):

$$\begin{array}{}\text{(5)}& {p{\mathrm{CO}}_{\mathrm{2}}}_{\mathrm{in}}={p{\mathrm{CO}}_{\mathrm{2}}}_{\mathrm{eq}}\exp \left[\mathrm{0.0423}\right(T_{\mathrm{in}}-T_{\mathrm{eq}}\left)\right],\end{array}$$

where ${p{\mathrm{CO}}_{\mathrm{2}}}_{\mathrm{in}}$ is the temperature-corrected partial pressure of CO₂, ${p{\mathrm{CO}}_{\mathrm{2}}}_{\mathrm{eq}}$ is the observed partial pressure of CO₂, T_in is seawater temperature measured close to the inlet and T_eq is the temperature measured in the tubing right before the equilibration chamber. Additionally, sea bottom temperature at 23 m depth below the inlet is measured by an internal thermometer inside an acoustic Doppler current profiler (ADCP).

The atmospheric concentration of CO₂ at the height of 57 m was measured at the atmospheric station run by Integrated Carbon Observation System, ICOS (C in Fig. 1). At this station, the dry molar fraction of CO₂ is measured using a very stable cavity ring-down spectroscopy technique (G2401, Picarro). For more information on the atmospheric ICOS measurements at Utö, the reader is referred to Kilkki et al. (2015).

3 Results and discussion

3.1 Suitability of the site for eddy covariance measurements

3.1.1 Environmental conditions

The measurement period (1 July–1 November 2017) represents the late summer and autumn seasons of the carbon cycle, when biological activity diminishes and the sea–air pCO₂ difference shifts from negative to positive (Fig. 3e). At the beginning of this period, from July to mid-August, the atmospheric pCO₂ exceeded that in the sea, and the sea acted as a sink of atmospheric CO₂. The average sea–air pCO₂ difference from July to mid-August was 8.0 Pa. During the latter part of August, the sea turned into a source as a result of diminished primary production. The efflux was enhanced when thermal stratification broke down and CO₂-rich water surfaced at the beginning of September, resulting in increased surface pCO₂. After this event, the partial pressure difference stayed predominantly within 10–20 Pa for the rest of the measurement period, showing only occasional deviations from this range.

The values reported for the sea surface pCO₂ represent conditions at the depth of 4.5 m at a single point 250 m away from the flux tower, within the flux measurement sector. In some cases, this location may not represent the carbonate conditions of the whole sea sector area. Especially during the beginning of the measurement period, it is possible that there are horizontal differences in the surface pCO₂ due to variations in biological activity (Rutgersson et al., 2008). During other periods, the horizontal differences are likely to be small in the open sea areas, and the surface sea layer is likely to be well mixed at least to the depth of 4.5 m. Thus, we assume that our pCO₂ observations represent the surface conditions.

Figure 3 Environmental conditions at Utö in July–November 2017: (a) wind speed, (b) wind direction (red band indicates the directions used for sea–air flux measurements in this study, and yellow band indicates additional directions that will be considered in the future), (c) stability (z∕L), (d) air temperature (red) and sea temperature at 4.5 m depth (orange) and at 23 m depth (blue), and (e) CO₂ partial pressure in air (red) and CO₂ partial pressure in seawater at 4.5 m depth (blue). Wind speed, wind direction, stability and air temperature were measured in the flux station (A in Fig. 1). Seawater temperatures and pCO₂ were measured at the inlet (B in Fig. 1). The pCO₂ in air was measured at the ICOS station (C in Fig. 1).

During our study period, the wind blew from the sea sector (180–260^∘) for 28 % of the time (Fig. 3). This caused gaps in the sea–air flux time series; for instance during the last 2 weeks of September, the wind directions were unsuitable for sea–air exchange measurements. The average wind speed within the marine wind sector was 6.9 m s⁻¹. Stability was mostly near neutral. The vertical rotation angle within the sea sector during the measurement period varied within 0–5^∘ with a mean of 3.1^∘ and a standard deviation of 1.4^∘, indicating that the flow divergence due to the cliff is limited. A small part (3.7 %) of the momentum fluxes measured in the marine sector were positive, and thus the corresponding CO₂ flux data were discarded.

In July, air temperature was still gradually increasing, and it reached 18 ^∘C by the end of the month. From the beginning of August, both air and sea surface temperature decreased approximately linearly. During the last 3 months, surface seawater cooled by 7 ^∘C and air by 10 ^∘C.

Figure 4 Normalized cospectral densities ${\text{Co}}_{w{\mathrm{CO}}_{\mathrm{2}}}$ of (a) the test setup and (b) the standard setup, and (c) Co_wT (green) and Co_wq (yellow). The solid lines represent the median cospectra, and the dashed lines represent the 10th and 90th percentiles. The black line is the model Co_wT (Horst, 1997). (d) Ratios of cospectral densities (dots) and the fitted exponential transfer functions (solid lines), where red indicates the standard setup, and blue indicates the test setup.

3.1.2 Spectral characteristics

The cospectra ${\text{Co}}_{w{\mathrm{CO}}_{\mathrm{2}}}$ of both measurement setups agreed well with the modeled cospectrum (Horst, 1997) in all frequencies (Fig. 4). The low-frequency ends of the measured cospectra approach zero in a similar way to the modeled cospectrum.

The attenuation of the highest frequencies was slightly higher in the more complex tubing: ${\text{Co}}_{w{\mathrm{CO}}_{\mathrm{2}}}$ of the test setup (with a drier and a virtual impactor) diverged from the modeled Co_wT at a slightly lower frequency than that of the standard setup. The half-power frequency for the standard and test system was 0.78 and 0.61 Hz, respectively. These values were used to correct the CO₂ fluxes according to the model cospectrum. The average high-frequency correction of the CO₂ flux during the measurement period was 16 % for the test setup and 12 % for the standard setup.

The drier and virtual impactor added approximately 1 m to the length of tubing, which only has a minor effect on the attenuation of turbulent fluctuations. However, the virtual impactor forms a 90 ^∘ angle to the tubing, which can stabilize the flow, making it laminar at a higher Reynolds number, Re, than in a straight tube (Lenschow and Raupach, 1991). It would be ideal to have turbulent conditions in the tubes, because scalar fluctuations dampen less in a turbulent pipe flow due to the more uniform velocity profile (Lenschow and Raupach, 1991). At low temperatures (0 ^∘C), Re of the standard setup was calculated to be 2980, whereas in the test setup it was 2380. By using the tube attenuation equations of turbulent flow presented by Massman and Ibrom (2008), the half-power frequencies were 13.7 and 7.2 Hz for the standard and test setup, respectively. Thus, at low temperatures, fluctuation attenuation due to transport in the tubes is not likely to affect the measurements. However, at higher temperatures (20 ^∘C), the Reynolds number of the test setup (2100) falls into the laminar flow region, because kinematic viscosity increases with increasing temperature. By using the tube attenuation equations for laminar flow by Lenschow and Raupach (1991), the half-power frequency of the test setup at 20 ^∘C was only 1.4 Hz, which is still higher than the empirically determined f₀.

The high-frequency attenuation of H₂O signals of the standard setup was corrected in order to accurately compare the CO₂ flux measurements with different setups as a function of latent heat flux. The high-frequency attenuation of H₂O flux of the standard setup was marked: f₀ ranged between 0.04 and 0.22 Hz, decreasing with increasing RH, which varied within 48–100 %. This attenuation is caused by sorption and desorption of water vapor at the walls of the tubings (Massman and Ibrom, 2008). During typical humidity conditions (RH=80 %), the flux attenuation is 44 % of the real H₂O flux. The ${\text{Co}}_{w{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ of the test setup was found to be distorted in all frequencies.

3.1.3 Stationarity

In order to evaluate the quality of our measurements, several theoretical assumptions of the eddy covariance method were analyzed. The flux footprint area was found to be horizontally homogeneous (Appendix B), and turbulence was well developed as integral turbulence characteristics could be expressed as functions of stability parameter (Appendix C). In this section, we analyze the fulfillment of the stationarity assumption and determine the effect of CO₂ flux nonstationarity on the CO₂ fluxes.

Over half of the flow conditions for which the CO₂ flux was determined, 65 %, were found to be nonstationary; i.e., the difference between the 30 min flux and the mean of the corresponding 5 min fluxes exceeded 30 %. Overall, 91 % of data were discarded due to unsuitable wind direction, nonstationarity or a positive momentum flux.

Figure 5 The effect of removal of nonstationarity cases on the CO₂ flux data: (a) no stationarity limit, (b) ${\text{RN}}_{\stackrel{\mathrm{‾}}{w^{\prime }{c}^{\prime }}}<\mathrm{0.6}$ and (c) ${\text{RN}}_{\stackrel{\mathrm{‾}}{w^{\prime }{c}^{\prime }}}<\mathrm{0.3}$. Red circle refers to the standard setup, and blue dot to the test setup.

We found that this stationarity criterion removed unphysical values effectively (Fig. 5). During July, the partial pressure difference was negative (−9.3 Pa on average), and thus downward (negative) fluxes were expected. Relaxing the stationary limit from 30 % to 60 % had only a small effect on data filtering. When no stationarity limit was applied, 20 % and 25 % of CO₂ sea–air fluxes that were measured during July 2017 by the test and standard setup, respectively, and were expected to be negative were positive. When the 60 % stationarity limit was used, only 5 % (test) and 6 % (standard) of the measured CO₂ fluxes were positive. With the 30 % stationarity limit, these numbers were only slightly lower: 3 % and 4 %, respectively.

Similarly high rejection rates have been obtained in previous studies. Miller et al. (2010) reported that 65 % of the sea–air CO₂ flux measurements failed a stationarity test in which the 13.7 min fluxes were compared with the averages of 3.4 min fluxes, and they assumed that CO₂ concentration heterogeneity could be a reason for the high rejection rate. Likewise, Blomquist et al. (2014) concluded that horizontal CO₂ concentration gradients caused by continental pollution sources can reduce the number of stationary situations. At Utö, concentration gradients may be produced by the horizontal heterogeneity due to the land–sea interface and the variations in marine primary production.

We examined how the absolute change in CO₂ concentration during the averaging period of 30 min relates to nonstationarity of CO₂ fluxes (Fig. 6). The change in CO₂ concentration was calculated as a difference in the mean CO₂ molar fraction of the last 30 s and the first 30 s. Typically, a change larger than 1 ppm was associated with rejection of the measurement due to the nonstationarity of CO₂ flux. Only 14 % and 18 % of the sea–air CO₂ fluxes for the test and standard setup, respectively, passed the stationarity test with the 30 % limit, for the cases with a CO₂ concentration change larger than 1 ppm. For the cases where CO₂ concentration changed more than 2 ppm during the averaging period, only 4 % and 6 % of sea–air CO₂ fluxes passed this stationarity test.

Figure 6 Relative nonstationarity of CO₂ flux as a function of the absolute change in CO₂ concentration during the averaging period: standard (red cross) and test (blue circle) setups. The green dashed line indicates the limit of ${\text{RN}}_{\stackrel{\mathrm{‾}}{w^{\prime }{c}^{\prime }}}=\mathrm{0.3}$.

If the timescale of the processes generating nonstationarity were less than the averaging time, it would be possible to reduce the amount of discarded nonstationary data by shortening the averaging time. To test for nonstationarity with shorter averaging periods, we calculated CO₂ fluxes for 15 min periods, which were compared with the average of 2.5 min CO₂ fluxes. It was found that this increased the number of accepted data only a little: by 6.1 % for the standard setup and 1.9 % for the test setup.

3.2 Setup comparison

The observed direction of the sea–air CO₂ fluxes was mainly consistent with the pCO₂ difference between the sea and the atmosphere, and both setups showed similar fluxes most of the time (Fig. 7). The standard setup tended to show slightly more positive fluxes than the test setup. In July, the pCO₂ difference was negative, and the mean sea–air CO₂ flux measured with the standard and test setup was −0.163 and −0.226 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, respectively. During the latter part of August (15–31 August), when the pCO₂ difference was shifting from negative to positive, the mean flux measured with the standard system (−0.050 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) was close to zero, while the test system still showed clearly negative fluxes (mean −0.122 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$). During August, the highest latent heat fluxes of the measurement period were measured, peaking at 120 W m⁻². During September, low latent heat fluxes (3.3 W m⁻² on average) were observed and both measurements showed closely similar monthly sea–air CO₂ flux: 0.188 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ (test) and 0.158 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ (standard). The highest absolute monthly CO₂ fluxes were measured in October, when pCO₂ difference was continuously high (16.6 Pa on average) and wind speed peaked occasionally, resulting in monthly sea–air CO₂ flux averages of 0.524 and 0.434 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ by the standard and test setup, respectively.

Figure 7 Sea–air fluxes of (a) CO₂ and (b) water vapor: standard setup (red circles) and test setup (blue crosses). Only the fluxes originating from the southwestern open sea sector (Fig. 1) are shown.

Figure 8 Measured gas transfer velocities as a function of wind speed. The orange line is the parametrization by Wanninkhof (1992).

As no previous measurement data are available for comparison, we validated the magnitude of the measured air–sea CO₂ fluxes roughly by calculating the gas transfer velocity ($k_{\mathrm{660}}=\frac{F}{K_{\mathrm{0}}\mathrm{\Delta }p{\mathrm{CO}}_{\mathrm{2}}}\sqrt{\frac{Sc}{\mathrm{660}}}$) from our data and compared it with the universal parametrization proposed by Wanninkhof (1992). Schmidt number (Sc) and solubility (K₀) of CO₂ were calculated according to Wanninkhof (1992) and Weiss (1974), respectively. For this comparison, we only included cases in which the absolute partial pressure difference (ΔpCO₂) between the sea and atmosphere was larger than 3 Pa. The gas transfer velocities derived from our measurements were mostly in good accordance with predictions (Fig. 8), which lends credence to our flux measurements. During medium wind speeds (3–10 m s⁻¹), the measured k₆₆₀ values line up closely with the parametrization. More measurements during low (<3 m s⁻¹) and high (>10 m s⁻¹) wind speeds are required to evaluate the relationship between winds peed and k₆₆₀ during these wind speeds.

Figure 9 The sea–air CO₂ flux measured with the standard setup vs. the test setup (equipped with a drier and a virtual impactor). The marker color indicates the concurrent latent heat flux, and the dashed line represents the 1:1 relationship.

Figure 10 CO₂ flux difference (standard setup–test setup) as a function of latent heat flux for (a) positive and (b) negative water vapor fluxes. The lines show a fitted mean together with the 10th and 90th percentiles for eight logarithmically equally spaced bins.

The Nafion drier eliminated the water vapor fluctuations effectively. The latent heat flux (measured with the standard system) was mostly positive, ranging within −28 to 122 W m⁻², with an average of 27 W m⁻². The H₂O variance measured with the standard setup was 0.021 mmol² mol⁻² on average, whereas the test setup measured the mean H₂O variance of 0.001 mmol² mol⁻². Even though the Nafion drier in our setup did not remove water vapor completely, it attenuated the fluctuations at all frequencies, in the same way that the tubing of a closed-path system attenuates temperature fluctuations. Water vapor molar fraction measured by the test setup varied within 2.9–8.8 mmol mol⁻¹, whereas the standard setup showed water vapor values of 5.7–15.7 mmol mol⁻¹. Thus, a total removal of water vapor from the sample air is not required to eliminate the water vapor fluctuations.

A high correlation was found between the two measurement setups (Fig. 9). The Pearson product–moment correlation coefficient between the sea–air CO₂ fluxes measured by the setups was 0.97. During low latent heat fluxes (<30 W m⁻²), the correlation coefficient increased to 0.99. As latent heat flux increases, the standard setup shows more positive CO₂ fluxes than the test setup (Fig. 10). These results are in agreement with the conclusions of Landwehr et al. (2014), who found that during low latent heat fluxes the difference between the dried and undried sea–air CO₂ flux measurements was very small. This also proves that a Nafion drier (in our case combined with a virtual impactor) does not disturb the CO₂ flux measurements.

A comparison of the CO₂ flux difference as a function of latent heat flux shows that the error is negligible for very small latent heat fluxes (<10 W m⁻²). With higher latent heat fluxes (∼60 W m⁻²), the mean difference is approximately 0.125 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$. Since the sea–air fluxes of CO₂ are typically small, the effect of water vapor on CO₂ flux measurement can cause a change in its sign, as observed during the latter part of August. However, as the mean difference increases, the scatter also increases. There are only few data points with negative latent heat fluxes, so no conclusions can be drawn concerning the effect of negative water vapor fluxes on CO₂ fluxes.

We showed that the use of the combination of a virtual impactor and a Nafion drier did not disturb the measurement of CO₂ fluxes. Closed-path instruments are typically protected by one or two Teflon filters (one close to the inlet and another next to the instrument), which should prevent any liquid water from reaching the instrument. However, the use of a Teflon filter close to the inlet may be unpractical, as a regular change of the filter in high-sea-spray conditions may prove laborious. In such a case, we suggest a virtual impactor as an option for the protection of the instrument and the tubing from water and sea salt.

4 Conclusions

In this paper, we presented new sea–air CO₂ flux measurements, comparing two closed-path gas analyzer setups installed on a shore of an island in the Baltic Sea. One of the setups was equipped with a drier and a virtual impactor. By inspecting several theoretical assumptions of the eddy covariance method, we showed that this land-based flux site is capable of effectively monitoring CO₂ fluxes between the atmosphere and the sea. The w–CO₂ cospectral densities calculated from our data showed an expected behavior with a moderate high-frequency loss. Turbulence was well developed, and the integral turbulence characteristics followed unique functions of stability and friction velocity.

We found that the two closed-path infrared gas analyzer setups showed similar CO₂ fluxes when latent heat fluxes were small. The correlation coefficient increased from 0.97 to 0.99 when only the CO₂ fluxes during small latent heat fluxes (<10 W m⁻²) were included. A higher latent heat flux resulted in a positive bias to the undried measurement. During high latent heat fluxes (∼60 W m⁻²), the difference between the CO₂ flux measurements was 0.125 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ on average, which is comparable to the monthly mean sea–air CO₂ flux. In July, the CO₂ partial pressure in the atmosphere exceeded that in the surface seawater, resulting in negative sea–air CO₂ fluxes, with monthly averages of −0.16 and −0.23 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ for the standard and test (drier) setup, respectively. In October, when the surface seawater had a higher CO₂ partial pressure than the atmosphere, the average sea–air CO₂ fluxes were 0.52 and 0.43 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ for the standard setup and test setup, respectively.

Our analysis indicated that, to obtain high-quality sea–air CO₂ flux data, a large proportion (65 %) of the original observations had to be discarded for nonstationarity. If the nonstationary cases were not discarded, 20 %–25 % of the sea–air CO₂ fluxes in July had a wrong sign, whereas only 3 %–4 % of the measured fluxes had an unphysical direction if the nonstationary cases were rejected. This nonstationarity was found to be linked to the changes in atmospheric CO₂ concentration during the averaging period. A change larger than 2 ppm was associated with a rejection rate of 94 %–95 % due to nonstationarity. The occurrence of nonstationary conditions was random in time, whereas unsuitable wind directions produced long, continuous gaps in the measurement time series. For the calculation of a CO₂ budget, these gaps must be filled, for instance by using the measured pCO₂ difference and a gas transfer velocity parameterized as a function of wind speed. Thus, continuous measurements of surface seawater and atmospheric CO₂ concentration at the Utö station provide useful additional data for the calculation of sea–air CO₂ flux exchange.

We showed that the use of the combination of a virtual impactor and a Nafion drier did not disturb the CO₂ flux measurement. While this configuration generated a slightly higher flux loss, we opt for this alternative for its reduced water vapor cross-sensitivity and better protection against sea spray.

Appendix A: Virtual impactor

The virtual impactor is based on two perpendicular air streams, which provide a means to separate particles by size (Fig. A1). The sample air stream is divided into minor and major flows: smaller particles are diverted to the major flow, and large particles with higher inertia, in this case water droplets, continue to the minor flow. The 50 % cut point of our virtual impactor was calculated to be 1.2–1.3 µm, i.e., a particle of this size has a 50 % probability of removal, while smaller particles are more likely to enter the major flow.

Figure A1 Schematics of the virtual impactor attached to the test setup. OD stands for outside diameter. The red pipes are enclosed by the black outer case. A more detailed figure is available on request.

Appendix B: Horizontal homogeneity

Horizontal homogeneity within the sea sector (180–340^∘) was examined in terms of the surface roughness length (z₀), which was calculated from the logarithmic wind profile law using data collected during neutral conditions, $|z/L|<\mathrm{0.01}$ (Fig. B1).

Figure B1 Surface roughness length as a function of wind direction. The red band indicates the directions that are considered the open sea sector used in this paper, and the yellow band indicates a wind sector that potentially provides useful data that were excluded from the present study.

The sea sector is clearly visible in the expected range of 180–340^∘, where z₀ is mainly lower than 1 mm. Outside this sector, z₀ is clearly higher, up to 1 m in northeastern wind directions. The sea sector indicated by z₀ coincides well with the sector deduced from the geographical map (Fig. 1).

While the sea sector can be considered sufficiently homogeneous, there is scatter in z₀ reflecting the fact that it depends on the sea state. Taylor and Yelland (2001) found that over a sea surface z₀ depends on the significant wave height and the slope of the wave. Since the measurement site is located on the coast, both the fetch and swell can depend on wind direction.

Appendix C: Turbulence

The Monin–Obukhov similarity theory predicts that statistical turbulence variables within the atmospheric surface layer are unique functions of the stability parameter (z∕L), which combines information about the height above ground, surface shear stress, surface heat flux and buoyant processes. These so-called integral turbulence characteristics have been found to have the same stability dependence over the sea as on land (Smith and Anderson, 1984) and can be used to examine the development of turbulence.

Figure C1 Integral turbulence characteristics. (a) ${\mathit{\sigma }}_w/u_{*}$ as a function of stability parameter z∕L, (b) ${\mathit{\sigma }}_w/u_{*}$ as a function of $z{f}_c/u_{*}$, (c) ${\mathit{\sigma }}_T/T_{*}$ as a function of z∕L during unstable stability and (d) ${\mathit{\sigma }}_T/T_{*}$ as a function of z∕L during stable stability.

Our data are in accordance with the results of Panofsky et al. (1977), who determined the relationship between the normalized standard deviation of vertical velocity and stability parameter:

$$\begin{array}{}\text{(C1)}& {\displaystyle \frac{{\mathit{\sigma }}_w}{u_{*}}}=\mathit{\alpha }\cdot (\mathrm{1}-\mathit{\beta }\cdot z/L{)}^{\mathrm{1}/\mathrm{3}},\end{array}$$

where u_* is friction velocity, α=1.3 and β=3.0. With large $-z/L$ values, the power law of 1∕3 fits our observations well, and the observed $\frac{{\mathit{\sigma }}_w}{u_{*}}$ approaches a constant of 1.3 in the neutral range similarly to Eq. (C1). Most of the observations fall into the stability range between −0.1 and −0.01. Within this range, the observed $\frac{{\mathit{\sigma }}_w}{u_{*}}$ is mainly scattered between values 1.0 and 1.5.

However, Monin–Obukhov similarity theory does not consider all the necessary information for describing turbulence characteristics in all conditions; e.g., the effect of Coriolis and pressure-gradient forces are excluded. Högström (1990) showed that during neutral stratification $\frac{{\mathit{\sigma }}_w}{u_{*}}$ can be described with z; u_*; and the Coriolis parameter, f_c. Smedman (1990) found that this relationship is valid at widely differing sites, including both terrestrial and marine surfaces. Thus we tested our data with neutral stratification, $|z/L|<\mathrm{0.1}$, against the function proposed by Högström (1990):

$$\begin{array}{}\text{(C2)}& {\displaystyle \frac{{\mathit{\sigma }}_w}{u_{*}}}=\mathit{\gamma }\cdot \log (z{f}_{\mathrm{c}}/u_{*})+\mathit{\delta },\end{array}$$

where γ=0.11 and δ=1.93. This function provides a reasonable fit to our data and a better description of $\frac{{\mathit{\sigma }}_w}{u_{*}}$ in the stability range between −0.1 and −0.01 than Eq. (C1) (Fig. C1).

For temperature, we compared the observations to the model of Foken et al. (1991):

$$\begin{array}{}\text{(C3)}& {\displaystyle \frac{{\mathit{\sigma }}_T}{T_{*}}}=\left\{\begin{array}{ll}\mathrm{1.0}\cdot {\left(\left|\frac{z}{L}\right|\right)}^{-\mathrm{1}/\mathrm{3}}& ,\frac{z}{L}\le -\mathrm{1}\\ \mathrm{1.0}\cdot {\left(\left|\frac{z}{L}\right|\right)}^{-\mathrm{1}/\mathrm{4}}& ,-\mathrm{1}\le \frac{z}{L}\le -\mathrm{0.0625}\\ \mathrm{0.5}\cdot {\left(\left|\frac{z}{L}\right|\right)}^{-\mathrm{1}/\mathrm{2}}& ,-\mathrm{0.0625}\le \frac{z}{L}\le \mathrm{0.02}\\ \mathrm{1.4}\cdot {\left(\left|\frac{z}{L}\right|\right)}^{-\mathrm{1}/\mathrm{4}}& ,\frac{z}{L}\ge \mathrm{0.02}\end{array}\right..\end{array}$$

Our observations were in accordance with the shape of this relationship, except for the near-neutral stability range, where the observations differ greatly from the model (Fig. C1c–d). As these cases have a very small sensible heat flux, $\frac{{\mathit{\sigma }}_T}{T_{*}}$ could be biased by division by a small number. However, this is not observed as the slope is well organized, and the spread is small. Our data suggest a steeper relationship should be used in this near-neutral stability range ($|z/L|<\mathrm{0.02}$):

$$\begin{array}{}\text{(C4)}& {\displaystyle \frac{{\mathit{\sigma }}_T}{T_{*}}}=\mathrm{0.1}\cdot {\left(\left|{\displaystyle \frac{z}{L}}\right|\right)}^{-\mathrm{0.94}}.\end{array}$$

On a logarithmic scale, the slope (−0.94) in our fit is almost twice as steep as the corresponding slope determined by Foken et al. (1991), i.e., −0.5.

Overall, our observations were well organized as a function of stability parameter, or friction velocity in the neutral case, and thus we did not discard any of the data based on the integral turbulence characteristics.