Heterogeneity Matters: Aggregation Bias of Gas Transfer Velocity Versus Energy Dissipation Rate Relations in Streams.

The gas transfer velocity, k , modulates gas fluxes across air-water interfaces in rivers. While the theory postulates a local scaling law between k and the turbulent kinetic energy dissipation rate ε , empirical studies usually interpret this relation at the reach-scale. Here, we investigate how local k ( ε ) laws can be integrated along heterogeneous reaches exploiting a simple hydrodynamic model, which links stage and velocity to the local slope. The model is used to quantify the relative difference between the gas transfer velocity of a heterogeneous stream and that of an equivalent homogeneous system. We show that this aggregation bias depends on the exponent of the local scaling law, b , and internal slope variations. In high-energy streams, where b > 1 , spatial heterogeneity of ε significantly enhances reach-scale values of k as compared to homogeneous settings. We conclude that small-scale hydro-morphological traits bear a profound impact on gas evasion from inland waters.

Introduction

Freshwater systems are a key component of the global carbon and nitrogen cycle (Battin et al., 2009; Bernhardt et al., 2018; Galloway et al., 2004; Hotchkiss et al., 2015; Marx et al., 2017; Raymond et al., 2013). In particular, headwater streams not only represent an important source of nutrients for downstream ecosystems but—owing to the high energy of the flow—they are also responsible for significant green‐house gas (e.g., , , and ) emissions to the atmosphere (Crawford et al., 2014; Kroeze et al., 2005; Hall & Ulseth, 2020; Horgby et al., 2019; Marzadri et al., 2017; Schelker et al., 2016). Thus, an accurate quantification of water‐air gas exchanges in low‐order streams proves crucial to better constrain current estimates of nutrient loads to recipient water bodies and river outgassing at regional or global scales (Rawitch et al., 2019).

One of the key parameters that control gas evasion in streams and rivers is the transfer velocity , which represents the water depth that equilibrates with the atmosphere per unit time. regulates the relation between gas fluxes across the air‐water interface and the concentration gradient across the water boundary layer, as per the Fick's first law of diffusion (Vingiani et al., 2021; Wanninkhof et al., 2009). The gas transfer velocity is a local property of the flow field, which reflects the intensity of the near‐surface turbulence and renewal rates in the water column (Zappa et al., 2003). Periodically, in fact, the fluid from the water column replaces the liquid film on the free surface where the exchange of gas with the atmosphere occurs. In streams and rivers, this circulation is typically related to the small eddies generated by the local turbulence at the Batchelor‐scale, in which is proportional to , where is the Schmidt number, is the kinematic fluid viscosity, is the molecular diffusivity of the dissolved gas, and is the energy dissipation rate per unit mass. In the above formulation, the Batchelor length scale, and thus gas evasion, depend on turbulence characteristics, gas diffusivity and temperature via and . The power‐law scaling of with is well established in the literature, and the theoretical value of the scaling exponent (0.25) is supported by many in‐field and laboratory experiments (Katul & Liu, 2017; Lamont & Scott, 1970; Moog & Jirka, 1999a). However, recent works (Ulseth et al., 2019) have proposed the existence of a super‐linear scaling between and in high‐energy streams, which arguably reflects bubble‐mediated processes promoted by air entrainment (Chanson et al., 2000).

Most of the available experimental methods to measure in the field—such as tracer experiments, night‐time regression of dissolved oxygen data, coupled estimation of daily metabolic rates—provide information about reach‐scale values of the gas transfer velocity. Consequently, available estimates of subsume spatial scales ranging from several tens to hundreds of meters—with shorter lengths that are typically associated with high‐energy streams, where gas exchange is enhanced (Bernot et al., 2010; Bott et al., 2006; Melching & Flores, 1999; Mulholland et al., 2001; Tsivoglou & Neal, 1976; Ulseth et al., 2019). While local estimates of are also feasible (e.g., via steady floating chambers, see Vingiani et al., 2021), reach‐scale estimates are particularly important for practical purposes, as long as they provide an objective basis for the upscaling of mass transfer rates at the network level (Bertuzzo et al., 2017; Raymond et al., 2013).

Owing to the spatial scales involved, experimental measurements usually refer to reaches that are featured by a strong internal heterogeneity. This lack of homogeneity of the flow field is particularly pronounced in high‐energy streams, which are often characterized by complex and diversified morphologies even at relatively short spatial scales (e.g., step and pools, slope changes, presence of hurdles within the flow, small cascades, dead‐end zones). Accordingly, experimental scaling laws between the mass transfer velocity and the energy dissipation rate available in the existing literature involve reach‐wise values of and —the latter being calculated from mean characteristics of the flow field (Moog & Jirka, 1999b; Raymond et al., 2012). Although empirical versus relations necessarily disregard the internal complexity of streams and rivers, the mechanistic link between the energy dissipation processes and the gas transfer velocity has an inherently local nature, and observed patterns of gas evasion within composite reaches properly reflect local heterogeneity of in‐stream hydrodynamic conditions (Kokic et al., 2018; Leibowitz et al., 2017; Vautier et al., 2020).

As the link between and is nonlinear (i.e., with ), spatial integrations along heterogeneous domains can alter significantly the nature of their relation (i.e., ), thereby making empirically derived laws potentially scale‐dependent. Accordingly, the main research hypothesis of the present work is that reach‐scale gas transfer versus total kinetic energy dissipation rates relations bear the signature of the spatial heterogeneity of the flow filed typical of many real‐world settings. This internal heterogeneity of river systems lies at the basis of a systematic aggregation bias affecting observed reach‐scale versus laws. This study combines a simple hydrodynamic model and a novel probabilistic framework to quantify and correct such bias using measurable external variables such as slope variations.

Theory

Hydrodynamic Setup and Local Turbulent Kinetic Energy Dissipation Rate

We focus on the one‐dimensional (1D) modeling of a rectangular cross‐sectional stream with a total length , here representing the length of an heterogeneous stream within which spatially integrated estimates of are available (e.g., tens to hundreds of meters—depending on the specific case‐study). Owing to the assumptions formulated hereafter, the minimum value of is ideally set to some tens of meters—although the general idea behind this study could be extended also to shorter reach segments. Regardless of the specific value of , we assume that owing to internal variations in the slope the turbulent kinetic energy dissipation rate (i.e., the power dissipated per unit mass by the flow) and the gas transfer velocity are continuous functions of —here representing the position along the considered stream path (Figure 1). To simplify the problem description, the impact of heat fluxes and temperature variations along the stream is here neglected, that is, . In addition, we also assume that the cross‐section width (), the discharge () and the bed roughness are constant along the focus reach. The hydraulic regime is then conceptualized as a sequence of piecewise nearly uniform flows, in which the possible transitions from supercritical to subcritical flow conditions (or vice versa) are neglected. In this setting, the changes in the bed slope induce local changes in the cross‐sectional velocity , according to the following relation: where is a positive exponent, is the water dept, and is a coefficient strictly related to the bed friction. In Equation 1, is the hydraulic radius, that is, the ratio between the cross‐sectional area and the wetted perimeter . For , simplifies as indicated in the r.h.s of Equation 1 because . The nature of the above equation sets the spatial resolution of the model to a minimum value of a few meters. As the discharge can be expressed as , we can rearrange Equation 1 and express the water depth, , as a function of the slope as: where is a coefficient strictly related to the particular empirical model chosen for Equation 1. For instance, if we specify Equation 1 according to the Manning's formula, we have (where is the Manning's coefficient) and , thereby leading to .

Figure 1

Schematic representation of a stream reach: map view (a), longitudinal profile of the elevation (b), and spatial patterns of energy dissipation rate, (Equation 3) (c). Panel (d) shows the probability density functions of within the reach (Equations 6 and 7).

According to Equations 1 and 2, both and vary with . This originates spatial variations in the energetic configuration of the flow, which are quantified by the turbulent kinetic energy dissipation rate of the bulk flow, . Under the assumption of uniform flow, can be approximated as (Tsivoglou & Neal, 1976): where is the gravitational acceleration. A word of caution is needed at this point. The proposed formulation neglects the vertical profile of turbulent kinetic energy dissipation rate within each cross section (Moog & Jirka, 1999b), and implicitly assumes that the value of driving the gas exchange at the water‐air interface can be approximated by the average value of in the bulk flow, as given by Equation 3. This choice is in line with the large majority of the experimental datasets currently available for streams and rivers (e.g., Raymond et al., 2012; Ulseth et al., 2019), and is suited to describe cases in which the gas exchange takes place over the entire water column, for example, owing to bubble‐mediated exchange and air entrainment.

Combining Equations (1), (2), (3), the following expression for is obtained: where is a constant, corresponding to the energy dissipation rate for a unit slope.

Characterizing Spatial Variations of Total Kinetic Energy Dissipation Rates

The spatial variability of is here described through its probability density function (pdf), which quantifies how often different values of are observed within the focus reach. In this framework, two different though related pdfs need to be introduced to quantify the spatial patterns of (Figure 1d):

The longitudinal pdf of , , which is the longitudinal pdf of the energy dissipation rate along the stream. This is simply the pdf of the values of observed within the different cross sections of the focus river reach. is the probability that belongs to the interval ()—which also corresponds to the fraction of length of the reach for which ; and

The volumetric pdf of , where is the bijective function that relates the water stage and the energy dissipation rate (see Equations 2 and 4) and is the mean stage along the focus reach. This is the volumetric pdf of the energy dissipation rate within the reach, and differently from attributes a larger weight to the values of that pertain to the regions of the domain where the stage is higher. Accordingly, is the fraction of water volume contained in the reach with a value of in the interval ().

Owing to the link between and (Equation 4), the two pdfs of introduced above can be calculated starting from the longitudinal pdf of the slope along the focus reach, . In line with observational evidence and previous studies, was here assumed to follow a two‐parameters Gamma distribution (Tarboton et al., 1989): where the shape parameter and the scale parameter constrain the moments of . In particular, is the inverse of the square of the coefficient of variation of the slope, that is, . The longitudinal pdf of the kinetic energy dissipation rates, can be obtained from Equations 4 and 5 as:

Likewise, the volumetric pdf of can be obtained combining the relation between and and Equation 6 as:

These pdfs express the probability distribution functions of the total kinetic energy dissipation rate within the focus reach, as a function of four parameters: two morphological parameters that define the internal heterogeneity of the slope ( and ), and two hydrodynamic parameters ( and ).

Upscaling Total Kinetic Energy Dissipation Rates

Let us consider the 1D stream reach of length depicted in Figure 1. The total power dissipated within the reach, , can be expressed by integrating along the reach the product between the dissipation rate , and the related water mass per unit length, that is,: where is the water density. The spatial patterns in the local energy dissipation rate originate a reach‐wise equivalent counterpart, which is here denoted as . The equivalent energy dissipation rate, , should provide the same total energy losses induced by the assortment of local experienced by the flow along the reach. The total dissipated power can be written in terms of as: where is the average cross sectional area of the stream reach. By equaling Equations 8 and 9, the equivalent energy dissipation rate can be written as: which corresponds to a weighted average of over the stream length, the weights being set by the volumetric pdf of in the considered reach, .

A more explicit formulation for can be obtained by inserting Equations 3 into 10 and recalling that , that is,: where is the mean water velocity measured by means of tracer injections.

Upscaling Gas Transfer Velocity

The concentration of a conservative gas in water is governed by the equation: where is the hydrodynamic dispersion, is the gas concentration in water, is the corresponding equilibrium concentration (the concentration that equilibrate with the atmosphere) and is the gas transfer velocity. If the dispersion term is neglected (i.e., ), as typically done in experimental tracer studies (e.g., Benson et al., 2014; Heilweil et al., 2016; Ulseth et al., 2019), the general solution of Equation 12 reads: where is the time‐dependent gas concentration in . Equation 13 shows how spatial variations of the local gas transfer velocity determine the patterns of downstream gas concentration. During field experiments, the internal heterogeneity of the studied stream is usually neglected and an equivalent reach‐scale gas transfer rate is estimated.

The relation between and can be found by equaling the concentration in the downstream section of the reach, , under the following conditions: (a) the case of an heterogeneous gas transfer velocity, , which corresponds to Equation 13 calculated for ; and (b) the case of a spatially uniform gas transfer velocity, which corresponds to Equation 13 calculated for with . This equality leads to the following expression for :

Equation 14 shows that the equivalent mass transfer velocity is a spatial longitudinal average of the local gas transfer velocities along the stream.

Upscaling Gas Transfer Velocity Versus Total Kinetic Energy Dissipation Rate Relations

Theoretical analysis and empirical results have long suggested the existence of a power‐law relation between and , of the type . While theoretical arguments suggest that , there is no consensus on the value of the scaling exponent and empirical data suggest that bubble‐mediated transport process could significantly increase up to 1.2 (Ulseth et al., 2019). Although these data mostly refer to reach‐scale estimates, they likely reflect the fact that the local scaling exponent might be higher than the theoretical value of 0.25 in case of air entrainment. For these reasons, here we treat as a model parameter. Regardless of the specific value of the scaling exponent, the relation has an inherent local nature implied by the mechanistic link among energy dissipation, turbulence, and gas exchange at the interface.

Nevertheless, such scaling laws have been empirically analyzed at the scale of entire reaches, including highly heterogeneous river segments. Using Equation 14, the equivalent mass transfer rate for a reach with length , , can be written as:

In Equation 15, the coefficient —which depends on the features of the underlying eddy dynamics—was assumed to be spatially uniform as we expect the turbulent structures to be weakly heterogeneous along the path. Equation 15 quantifies the equivalent gas transfer rate of an heterogeneous stream with a mean value of the total kinetic energy dissipation rate equal to . This equivalent mass transfer rate, , is equal to the order moment of , seen as spatial random variable with a longitudinal pdf . Were the local‐scale relation between and preserved at the reach level, the r.h.s. of Equation 15 should read . As in the general case , Equation 15 shows that the local scaling law between and is in fact altered when it is aggregated and re‐interpreted at the reach level, unless the pdf of is a Dirac delta distribution (i.e., the slope is spatially uniform within the reach). The resulting aggregation bias, is quantified as the relative difference between the reach scale gas transfer velocity of an internally heterogeneous stream (Equation 15) and the gas transfer velocity of an equivalent homogeneous stream with the same mean turbulent kinetic energy dissipation rate (). Under the assumptions made, the following analytical expression for can be obtained (see Supporting Information):

Interestingly, Equation 16 shows that the bias does not depend on the discharge, , and the friction coefficient, , but is only a function of the following independent parameters: (a) the exponent of the local scaling law that links the mass transfer rate to the total kinetic energy dissipation rate, ; (b) the hydraulic exponent , which quantifies the relation between the hydraulic radius and slope; and (c) the parameter , which quantifies the relative variations of the slope within the focus reach ().

In Figure 2, we have represented the aggregation bias as a function of the coefficient of variation of the slope and the exponent of the local scaling law in the case . Similar plots can be obtained for different values of . The plot shows that for the bias is positive, and increases nonlinearly for more pronounced internal variations of the bed slope. When the variability of the slope is limited (), is significant only when exceeds 1. When the internal variations of the slope are moderate (), instead, the aggregation bias is important only for . For values of in the range , exceeds 10% only when is above 1, whereas for the magnitude of the bias is usually negligible. Interestingly, for (which is the expected scaling exponent according to the theory), the bias is negative up to . The behavior shown in Figure 2 is generated by the combined action of two independent components that are responsible for generating the overall aggregation bias. The first component is termed volumetric bias component and provides a positive contribution to the bias due to the difference between longitudinal and volumetric averages of in a stream, with in nonhomogeneous reaches; the second component can be termed nonlinearity bias component and it relates to the difference between the order moment of a random variable and the power of its mean, with the latter being smaller (larger) than the former for (). Thus, a negative bias is observed only when and the bias is dominated by the nonlinearity component, which is negative in this case. In all the other instances, instead, the bias is positive, as shown by Figure 2.

Figure 2

Estimation bias as a function of the coefficient of variation and the exponent of the local scaling law, . The blue line refers to the theoretical scaling exponent of 0.25, the orange line refers to the observed scaling exponent in low energy streams as given by Ulseth et al. (2019) (0.35) and the green line refers to the observed scaling exponent in high‐energy streams as given by Ulseth et al. (2019) (1.18).

Overall, the analysis suggests that the aggregation bias can be neglected in low energy streams where the scaling exponent is small (theoretical value of , observed value of , see Zappa et al., 2007 and Ulseth et al., 2019). In those settings, the aggregation bias might be particularly low also because of the limited slope variations, which characterize most of the low‐energy streams where tracer gas experiments were performed, especially in cases when the analyzed reaches are relatively short. The limited effect of the aggregation bias for these systems is explained by the compensating effects associated with the two constitutive components of the bias. Instead, the aggregation bias is particularly important within high‐energy streams, where the scaling exponent is higher and the coefficient of variation of the slope can be of the order of 1. Therefore, a specific real‐world example pertaining to this class of streams will be further discussed with the aid of numerical simulations in the next section.

Numerical Simulations

A real‐world headwater catchment in the Italian Alps was analyzed to quantify the spatial patterns of , , , and the corresponding aggregation bias in a representative high‐energy stream. The selected case study is the Rio Valfredda, which drains a 5.3 km2 catchment in the Alps of northern Italy. The elevation ranges from 1,500 to 3,000 m a.s.l., generating significant spatial heterogeneity in lithology and soil cover. The alpine climate typical of the area is characterized by short rainy summers and long rigid winters with significant snowfall. The average annual precipitation is 1,500 mm unevenly distributed along the year, with low flows in winter and higher discharges during spring and summer. The stream network is 16 km long and reaches the fourth Strahler order. The channel slope averages at 26.7%, with maximum values of 125% in the northern and eastern parts of the network. The river bed composition is made up of rock emergencies, gravel, and silt, while the channel width ranges from 30 cm to 1.5 m. A LiDAR survey was carried out in 2018 to generate a Digital Terrain Model (DTM) with a 20 cm resolution. For more information about the catchment, the reader is referred to Durighetto et al. (2020) and Botter and Durighetto (2020).

The aggregation bias in the Valfredda creek was estimated by means of a numerical application of the hydrodynamic model within 200 reaches with a fixed length, , of 100 m and a random position along the network. The value of was chosen to represent the typical length of the reaches where tracer studies are performed in high‐energy settings, though alternative choices would have led to analogous results. For each reach, the longitudinal profile was reconstructed from the DTM, and the spatial patterns of local slope were calculated with a resolution of 2 m. Considering the bed composition of the Rio Valfredda (and the related roughness), the uniform flow is assumed to be supercritical within most of the analyzed stream network (Froude number is within 95 of the network length). Equations 2 and 4 were used to calculate and along each river reach, assuming . The local gas transfer rate was then calculated using the law proposed by Ulseth et al. (2019) for high‐energy streams: . Equations 10 and 14 were employed to calculate the reach‐wise values of and . Finally, the aggregation bias was calculated as .

Figure 3a shows as a function of the coefficient of variation of the slope, , for the 200 sample reaches of the Valfredda creek. The simulated aggregation bias follows the general trend predicted by Equation 16, with a nonlinear increase of with . In approximately 30% of the cases , with peaking at 80% when . The scatter of the numerical points around the theoretical equation is generated by deviations of the slope pdf from the gamma distribution. This happens more often in the most heterogeneous reaches, where the step and pool configuration generates a bimodal distribution of .

Figure 3

Estimation bias as a function of the coefficient of variation of the slope (a). Each point represents a different 100 m reach of the Valfredda stream network. Note that a maximum overlap among different reaches of 40% was allowed in this analysis. In this example, was set to 1.18 as suggested by Ulseth et al. (2019) for high‐energy streams ( is systematically higher than 0.02 m2 s−3). Spatial map of the slope of the stream network in the Valfredda catchment (b).

Discussion

The main goal of this study is to clarify analogies and differences between local and reach‐scale relations. To this aim, a simple hydrodynamic model was proposed in which the spatial variations of and are linked to the heterogeneity of the slope. While the formulation is quite general, the suitability of the proposed model to describe the water flow in mountain reaches needs to be carefully assessed case‐by‐case. In these settings, in fact, the flow can be tortuous and irregular, thereby challenging the use of a simplified 1D framework. Several studies suggested the use of an updated version of the classical Gauckler‐Manning equation to model water flows in small mountain streams in which the Manning coefficient varies with the hydraulic radius (see Marcus et al., 1992). The proposed model could incorporate the effect of this dependence of the roughness on by properly tuning the coefficients and in Equation 1. Furthermore, local energy dissipations with enhanced gas evasion could be observed in correspondence of hydraulic jumps (which are not described by our model), owing to the transition from supercritical and subcritical flow. While we acknowledge that and formulations might be more flexible in describing the spatial heterogeneity of complex flow fields typically observed in high‐energy streams (including those observed in presence of, for example, step and pools, hydraulic jumps, and abrupt planar discontinuities), we propose that the concept of aggregation bias introduced in this paper is quite general. In particular, our study indicates that—regardless of the specific features of the flow field—internal heterogeneity of the energy dissipation processes generates a potential bias in spatially integrated laws. The magnitude of the bias is eventually driven by the local scaling exponent and the underlying pdfs of the turbulent kinetic energy dissipation rate in the focus reach (Equations 10 and 14), which could be evaluated on a case‐by‐case basis. Therefore, the proposed framework not only provides a first‐order assessment of the impact of hydrodynamic heterogeneity on the outgassing of streams with variable slope, but offers a robust conceptual basis for evaluating the aggregation bias in many other settings, in which the major assumptions of the hydraulic model developed here are not fulfilled.

Our analytical results indicate that the aggregation of spatially heterogeneous flow conditions generates reach‐scale that are higher (or slightly lower) than the transfer rates that would be observed in a homogeneous system with the same mean turbulent kinetic energy dissipation rate if (or ). The magnitude of the aggregation bias, seen as the relative difference in the mass transfer rate of the homogeneous and the heterogeneous system for a given value of the mean turbulent kinetic energy dissipation rate, mainly depends on two factors: (a) the scaling exponent of the local versus law, with larger biases observed for ; and (b) the coefficient of variation of the local slope. In particular, for (which is the observed apparent scaling exponent for high‐energy streams) and (as frequently observed in steep headwaters) the bias exceeds 100%. These results are also confirmed by the numerical simulations performed within 200 reaches belonging to a small headwater catchment of the Italian Alps. The numerical results demonstrate the robustness of our analytical results and the reliability of the assumption made about the distribution of the slope. Overall, our analysis reveals that theoretical local relations cannot be straightforwardly inferred from aggregated versus data points in high‐energy streams (where ) unless the system is internally homogeneous, which is a very rare circumstance in most real‐world settings.

The practical implications of the theoretical results shown in the paper are manyfold. First of all, the apparent scaling exponent that might be deduced from multiple reach scale tracer studies in high‐energy streams could be different from the underlying scaling exponent of the local versus relation, especially whenever is a function of , that is, if more energetic streams are characterized by a different degree of internal heterogeneity as compared to less energetic streams. Therefore, regional estimates and available databases should be analyzed with caution as they could be affected by systematic aggregation biases. Second, reach‐scale estimates pertaining to streams characterized by a different degree of heterogeneity might not be directly comparable: in high‐energy streams, where empirical data suggest values of larger than 1, more heterogeneous streams will be featured by relatively larger gas transfer rates for a given value of the total kinetic energy dissipation rate. This might also be one of the reasons that underlie the observed scatter in the reach‐scale experimental versus pairs in Ulseth et al. (2019). Understanding how reach‐wise scaling laws originate from local relations is also important for spatial extrapolations of gas transfer velocities, which lie at the basis of regional or global‐scale studies on greenhouse gas evasion from streams and rivers. Our analysis indicates that spatial extrapolations of should be performed taking into account the internal heterogeneity of streams: river reaches with the same mean energetic level could be in fact characterized by highly heterogeneous gas transfer velocities, if they are internally different. Thus, small‐scale morphological heterogeneity of low‐order streams (e.g., steps, pools, and cascades) arguably represents an important missing component of in‐stream nutrient cycling, which has relevant implications for the quantification of gas evasion from inland freshwater systems.

Conclusion

In this study, a simplified hydrodynamic model was used to analyze how a local power‐law relation between mass transfer velocity and energy dissipation rate can be altered when it is spatially integrated along a heterogeneous stream, and then referred to equivalent reach‐scale values of and . Analytical derivations and numerical simulations indicate that local power‐law relations are not preserved at the reach (or network) level in presence of significant internal spatial variations of the energy dissipation rate. In particular, a systematic positive aggregation bias is observed when the exponent of the local versus relation, is larger than 0.5, as typically observed in high‐energy streams. In these settings, the presence of spatially heterogeneous flow conditions generates reach‐scale values of that are higher than the corresponding mass transfer rates that would be observed in a homogeneous system with the same mean turbulent kinetic energy dissipation rate. The magnitude of the bias can be as high as 100% if and the standard deviation of the slope exceeds the mean slope (i.e., ). On the other hand, our results indicate that the aggregation bias is negative and negligible when the scaling exponent of the versus relation is lower than 0.5, as predicted by surface renewal theories and empirically observed in low‐energy streams. Therefore, available datasets should be analyzed with caution as the presence of aggregation biases in high‐energy streams might impair their use for the quantification of global‐scale gas evasion rates from rivers.

Spatially averaged quantity.

Cross‐sectional area ().

Proportionality coefficient of the local relation ().

Power‐law scaling exponent of the local relation ().

Gas concentration in water ().

Gas concentration at the reach inlet ().

Equilibrium gas concentration in water ().

Molecular diffusivity ().

Hydrodynamic dispersion ().

Hydraulic coefficient related to the uniform flow equation ().

Relative aggregation bias ().

Gravitational acceleration ().

Local water depth (stage) ().

Gas transfer velocity ().

Reach‐wise equivalent gas transfer velocity ().

Reach length ().

Manning's coefficient ().

Total power dissipated along the reach ().

pdf of the reach slope ().

Longitudinal pdf of the turbulent energy dissipation rate ().

Volumetric pdf of the turbulent energy dissipation rate ().

Water discharge ().

Hydraulic radius ().

Local slope ().

Schmidt number ().

time ().

Local uniform flow velocity ().

Uniform reach width ().

longitudinal coordinate along the flow path ().

Scale parameter of the pdf of the slope ().

Shape parameter of the pdf of the slope ().

Local turbulent kinetic energy dissipation rate ().

Local turbulent kinetic energy dissipation rate for unit slope ().

Reach‐wise equivalent energy dissipation rate ().

Longitudinal average of ().

Volumetric average of ().

Gamma function.

Exponent of the hydraulic radius in the uniform flow equation ().

Kinematic viscosity ().

Friction coefficient ().

Water density ().

Wetted perimeter of the cross section ().

Supporting Information S1

Click here for additional data file.