Local processes with a global impact: unraveling the dynamics of gas evasion in a step-and-pool configuration

Headwater streams are important sources of greenhouse gases to the atmosphere. The magnitude of gas emissions originating from such streams, however, is modulated by the characteristic microtopography of the riverbed, which might promote the spatial heterogeneity of turbulence and air entrainment. In particular, recent studies have revealed that step-and-pool configurations, usually found in close sequences along mountain streams, are important hotspots of gas evasion. Yet, the mechanisms that drive gas transfer at the water–air interface in a step-and-pool configuration are not fully understood. Here, we numerically simulated the hydrodynamics of an artificial step-and-pool configuration to evaluate the contribution of turbulence and air entrainment to the total gas evasion induced by the falling jet. The simulation was validated using observed hydraulic features (stage, velocity) and was then utilized to determine the patterns of energy dissipation, turbulence-induced gas exchange and bubble-mediated transport. The results show that gas evasion is led by bubble entrainment and is mostly concentrated in a small and irregular region of a few square decimeters near the cascade, where the local gas transfer velocity (k) peaks at 500 m d⁻¹. The enhanced spatial heterogeneity of k in the pool does not allow one to define a priori the region of the domain where the outgassing takes place and makes the value of the spatial mean of k inevitably scale-dependent. Accordingly, we propose that the average mass transfer velocity should be used with caution to describe the outgassing in spatially heterogeneous flow fields, such as those encountered in step-and-pool rivers.

1 Introduction

Headwater streams play a pivotal role in the biogeochemical functioning of river networks and regulate key ecosystem services associated with fluvial environments (Appling et al., 2018; Bernhardt et al., 2018). Furthermore, mountain streams not only transport matter and chemicals from the uplands to downstream waterbodies but are also able to exchange gaseous compounds with the overlying atmosphere (Butman and Raymond, 2011; Battin et al., 2009). Despite the small flow rates and reduced water–air interfaces that characterize low-order streams, the efficiency of gas transfer in high-energy mountain settings makes these waterbodies a significant source of greenhouse gases (such as CO₂ and N₂O) to the atmosphere. Thus, enhancing our understanding of gas evasion from headwater streams would have important implications for a broad range of fields in biology and environmental sciences (Kroeze et al., 2005; Crawford et al., 2014; Schelker et al., 2016; Marzadri et al., 2017; Horgby et al., 2019; Hall and Ulseth, 2020).

Yet, quantifying the cycling and release of greenhouse gases in low-order streams remains challenging, owing to the intertwined effect of different biotic and abiotic agents. In particular, the physical mechanisms responsible for gas evasion in river networks are often modulated by the hydrodynamic features of the flow field, which are in turn strongly impacted by the (micro)topography of the riverbed (Duvert et al., 2018; Botter et al., 2021; Rocher-Ros et al., 2019; Looman et al., 2021). Differently from floodplain rivers – which exhibit smooth hydraulic conditions – the steep slopes and irregular beds typical of mountain streams yield complex and heterogeneous flow patterns. In this setting, the free surface of the water flow can frequently break, owing to obstacles, bends or abrupt variations in the channel bottom.

Recent studies evidenced that local steps are important hotspots of gas evasion, to the point that in all the settings where steps are distinct from the stream's morphology (e.g., in high-energy rivers or in streams originated by glacial processes), they control the overall network-scale outgassing (Vautier et al., 2020; Botter et al., 2022). While it has been argued that gas transfer associated with riverbed drops is promoted by the concurrence of high turbulence and bubble entrainment (see, e.g., Ulseth et al., 2019), there is limited knowledge of the physical processes that drive outgassing dynamics in local steps (see, e.g., Cirpka et al., 1993). Moreover, to date, the spatial patterns of energy dissipation and gas exchange in the plunging jet and the receiving pool of individual steps have not been analyzed. Consequently, the characteristic length scales of turbulence-induced and bubble-mediated gas transport in a step-and-pool configuration are still unknown. The characterization of the local and spatially heterogeneous nature of gas evasion in step-and-pool formations has important practical and theoretical implications. In particular, the spatial heterogeneity of gas exchange along a step constrains our ability to conceptualize such geomorphic elements by identifying a priori the size of the surface through which most of the gas exchange occurs (i.e., the step gas footprint) and the corresponding “effective gas transfer velocity” therein.

The present work aims at filling these knowledge gaps by combining experimental data and numerical modeling. In particular, here we reproduce the hydrodynamics of a simple step-and-pool geometry within an artificial flume using numerical simulations, which are then used to model the underlying gas exchange mechanisms taking place across the step. The numerical analysis allows us to (i) quantify the spatial variability in the gas transfer velocity, discerning the contributions of the turbulence and the entrained in-water bubbles, and (ii) identify the spatial patterns of energy dissipation and gas evasion in a step-and-pool configuration. The major implications of our findings for future gas transfer studies are also discussed.

2 Methods

2.1 In-flume experiment

An experiment to study the hydrodynamics of a step-and-pool configuration was carried out under controlled hydrodynamic conditions inside a horizontal, 6 m long plexiglass-made artificial flume with a rectangular section of width and height equal to 0.3 and 0.5 m, respectively. The water circulates through the channel via a constant head tank that maintains a steady discharge of 2 L s⁻¹, which is accurately measured by a magnetic flowmeter.

A step-and-pool formation was artificially created in the flume by inserting a broad-crested weir and a tapered positive step (see Fig. 1). The height of the two elements was 25 and 9 cm, respectively. Owing to the difference in the water level upstream and downstream of the weir, a local step with a water height drop ΔH equal to ≈ 15 cm was obtained. The pool, bounded by the weir and the step, had a depth of z₀≈ 11 cm. The pool was long enough (100 cm, about 10 times the water depth) so that the hydrodynamics near the jet were not impacted by the presence of the downstream positive step. Finally, an orifice on the wall of the weir ventilated the air below the cascade jet to maintain atmospheric pressure therein.

Figure 1 Schematic of a the step-and-pool configuration. (a) Lateral and (b) frontal view of the flume. Lengths are expressed in centimeters.

We measured the water level along the step-and-pool configuration using a graduated hydrometric rod (accuracy ± 0.1 mm) mounted on a moving carriage. Measurements were taken every 2 cm along the middle axis of the flume, allowing us to reconstruct the free-surface profile in the region of interest.

Local flow velocity was assessed in some representative sections along the middle axis of the pool by two acoustic Doppler velocimeters (ADVs). To this aim, we used a Vectrino Profiler ADV (Nortek, Rud, Norway) to record the velocity near the bottom. The profiler provides observations of the three-component velocity, u_i, with a resolution as fine as 1 mm over a 30 mm range placed 40 mm below the transmitter and an output rate of 100 Hz. For each section, we carried out a series of measures as follows: starting from the bottom, we repeatedly raised the ADV in 10 mm increments, recording velocities for 30 mm of water column each time. Since the ADV works correctly only when the transmitter is dipped in the water, the fluid volume close to the air–water interface (with a distance from the bed of more than 70 mm) was not analyzed. Therein, the velocity was assessed by a 2D SonTek ADV (Xylem Inc., Rye Brook, NY, USA). This ADV provides two-component velocity measures for a single point that corresponds to a control volume of 1 cm³, with a maximum sampling rate of 50 Hz. Measures of the horizontal components of the velocity were taken every 10 mm for 70 mm starting from the free surface.

To provide statistically consistent samples, the ADVs recorded the velocity signal for more than 30 s. During the tests, talcum powder was added to the water to achieve a signal-to-noise ratio (SNR) higher than 30 dB, thus minimizing the background noise recorded by the instruments.

2.2 Numerical simulation of the step-and-pool configuration

The in-flume experiment was numerically simulated using the software FLOW-3D Hydro 2022R1 (Flow Science Inc., Santa Fe, NM, USA), which adopts grids based on structured orthogonal cells to shape complex geometries using the fractional area–volume method (FAVOR). The evolution of the free surface is instead modeled using the volume-of-fluid (VOF) method (Hirt and Nichols, 1981), which can also track complex interfaces.

To simulate the hydrodynamic field of the step-and-pool configuration, the software solves the Reynolds-averaged Navier–Stokes equations along the Cartesian axes x_i coupled with the RNG (re-normalization group) e_k–ε scheme (e_k and ε being the turbulent kinetic energy per unit of mass and the turbulent kinetic energy dissipation rate, respectively), which resolves the turbulence at the sub-grid scale (Yakhot and Orszag, 1986).

The simulated dynamics in the pool included a full description of the bubble entrainment associated with the water flow. The air model used in this paper describes air transport in the water column according to the standard advection–diffusion equation of a scalar with known concentration C. In the simulation, the air bubbles were assumed to be spheres with a radius, r, equal to 0.5 mm, a value congruent with the analysis of cam-derived images taken during the experiment and in line with previous experimental studies (Woolf et al., 2007; Klaus et al., 2022; Karn et al., 2016). The model allowed us to quantify the entrained and escaping bubble fluxes across the free surface, q_B.

In the numerical experiment, the cell size ranged from 2 mm in the region where the jet plunged to 8 mm upstream of the broad crest weir. In total, the flume was discretized using approximately 7.5 million cells. The boundary conditions were the following: (i) a steady inflow of 2.0 L s⁻¹ upstream of the weir and (ii) free outflow downstream of the step. The no-slip velocity constraint was set on both the bottom and channel sides, where the equivalent roughness was set at ${\mathit{\delta }}_{\mathrm{e}}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{2}}$ mm. Initially, fluid was at rest, and transient conditions were observed for less than 20 s before stationary conditions were reached in the domain.

2.3 Modeling gas exchange in step-and-pool configurations

The gas exchange between water and the atmosphere takes place through the thin separation layer at the interface between the two fluids. The exchanged volumetric gas flux, F, is proportional to the difference in gas concentration between air and water and is described by Henry's law as

$$\begin{array}{}\text{(1)}& F=k\left(C_{\mathrm{w}}-{\displaystyle \frac{C_{\mathrm{a}}}{K_{\mathrm{H}}}}\right),\end{array}$$

where k is the gas transfer velocity (representing the water depth equilibrated with the atmosphere in the unit of time); C_w and C_a are the gas concentration in water and atmosphere, respectively; and K_H is the Henry constant determined for the specific exchange process under investigation. When both fluids are at rest, the exchange rate depends on the molecular diffusivity of the dissolved gas, D, which slowly drives the absorption to (or the release from) the water column. Instead, in turbulent flowing water k is much higher than that induced by the molecular diffusivity. In this case, the eddies originating from turbulence continuously renew the water close to the free surface, enhancing the transfer of gas at the interface (Lamont and Scott, 1970; Katul and Liu, 2017). Hereafter, the gas transfer velocity driven by turbulent mechanisms of this type is defined as k_T.

However, in many real-world settings, air entrainment and bubbles can intensify gas exchange in heterogeneous flow fields, primarily by increasing the total exchange area between water and air. The bubble-mediated gas exchange is encapsulated here by the bubble-mediated gas transfer velocity, k_B. Therefore, in turbulent flowing waters where air entertainment takes place, the total gas transfer velocity k can be computed as the sum of velocities due to free-surface and bubble-mediated exchange (Klaus et al., 2022):

$$\begin{array}{}\text{(2)}& k=k_{\mathrm{T}}+k_{\mathrm{B}}.\end{array}$$

The separation of the contributions of k_T and k_B to the total gas exchange in a step-and-pool configuration is one of the main contributions of this paper. The following subsections describe how these terms are computed in our numerical simulation.

2.3.1 Turbulence-induced air–water gas exchange in the flume

Existing estimates of gas exchange in the presence of turbulent flows are based on the idea of linking k_T with some key characteristic hydrodynamic quantity (e.g., the energy dissipation). In mountain streams, with low water depths and high velocities, the theoretical schemes that are better suited to quantify k_T are those based on small-eddy models (Moog and Jirka, 1999b). In this framework, the renewal of the free surface of the water column is led by the small-scale vortexes originated by the turbulence near the interface, which in turn are controlled by the rate of dissipation per mass unit of the turbulent kinetic energy, ε. According to this model, k_T scales with ε as follows:

$$\begin{array}{}\text{(3)}& k_{\mathrm{T}}={\mathit{\alpha }}_{\mathrm{T}}{\mathit{Sc}}^{-\mathrm{0.5}}{\left(\mathit{\nu }\mathit{\epsilon }\right)}^{\mathrm{0.25}},\end{array}$$

where α_T=0.2–0.4 is a calibration factor, and $\mathit{\nu }\approx \mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹ is the kinematic viscosity of water. Moreover, in Eq. (3), $\mathit{Sc}=\mathit{\nu }/D$ is the Schmidt number, which expresses the link between kinematic and molecular diffusivity. The nature of the processes described by Eqs. (1) and (3) is inherently local, and they occur at the Batchelor scale (Batchelor, 1959), which is defined as ${\mathit{\lambda }}_{\mathrm{B}}=\mathit{\eta }/\sqrt{\mathit{Sc}}$, where $\mathit{\eta }=({\mathit{\nu }}^{\mathrm{3}}/\mathit{\epsilon }{)}^{\mathrm{1}/\mathrm{4}}$ is the Kolmogorov length, i.e., the characteristic length of the microvortexes. In running waters, Sc>1, the spatial scale of the gas transport is thus smaller than that of the smallest eddies generated by turbulence in its energy cascade, that is, the process of the energy transfer from the large-scale eddies to microvortexes (Tennekes and Lumley, 1972).

Previous experimental studies have demonstrated that the exponent −0.5 in Eq. (3) applies to cases in which the free surface is smooth, while the same exponent may decrease to −0.67 if the free surface is riffled (Zappa et al., 2007; Jähne et al., 1987). It should be noted that, in the general case, the gas exchange also depends on the fluid temperature, which affects both ν and D. However, since the present work aims at reproducing a short portion of a stream in which the water temperature is nearly constant, we neglect the thermal effect on D and ν so that Sc is assumed to be constant, and the gas transfer velocity k_T scales with ε^0.25.

Nevertheless, in-field estimations of gas transfer velocity are usually performed at the reach scale based on spatially averaged quantities under the assumption of uniform flow. Under the above circumstances, the mean turbulent kinetic energy dissipation rate $\stackrel{\mathrm{‾}}{\mathit{\epsilon }}$ (hereinafter the overline applies to spatially averaged quantities) can be expressed as (Ulseth et al., 2019; Hall and Ulseth, 2020; Hall and Madinger, 2018; Raymond et al., 2012)

$$\begin{array}{}\text{(4)}& \stackrel{\mathrm{‾}}{\mathit{\epsilon }}=g\stackrel{\mathrm{‾}}{U}\stackrel{\mathrm{‾}}{S},\end{array}$$

with g the gravity acceleration, $\stackrel{\mathrm{‾}}{U}$ the cross-sectional flow velocity and $\stackrel{\mathrm{‾}}{S}$ the average bed slope. It should be noted that the proposed scheme is suitable if the turbulence generated at the bottom is the main driver of the overall intensity of the turbulence in the water flow. In headwater streams, instead, riverbed heterogeneity, partially emergent boulders and small cascades can be important additional turbulence sources (Moog and Jirka, 1999a; Botter et al., 2022). In the particular setup analyzed in this study, the plunging jet of a cascade dissipates its energy by inducing high ε in the outermost fluid layer, thereby preventing the use of Eq. (4).

Consequently, k_T is calculated based on the computed value of ε near the free surface, using Eq. (3) with α_T=0.4 and standardizing the result with a Schmidt number equal to 600.

2.3.2 Bubble-mediated gas exchange

In this work, while dealing with bubble-mediated gas exchange we focus on the so-called “kinematic bubble effect” (Liang et al., 2013). According to our approach, the amount of gas transferred from the water to the bubbles (or vice versa) is estimated using the following three-step procedure: first we determine the reference bubble lifetime, T_B (i.e., the time spent by the bubble within the water column when it reaches the surface; Step A); then, we calculate the gas transferred across the bubble interface according to the difference in concentration compared to the surrounding bulk flow along the bubble flow path (Step B); finally, the total gas transfer velocity is estimated by considering the total amount of bubbles involved in the process (Step C).

• Step A. A reasonable estimation of the bubble lifetime, T_B, is given by

  $$\begin{array}{}\text{(5)}& T_{\mathrm{B}}={\mathit{\alpha }}_{\mathrm{B}}{\displaystyle \frac{z_{\mathrm{0}}}{u_{\mathrm{B}}}},\end{array}$$

  where α_B is an O(1) coefficient, and u_B is the bubble rise velocity. Here, u_B was estimated using the model proposed by Woolf (1993):

  $$\begin{array}{}\text{(6)}& u_{\mathrm{B}}=\left\{\begin{array}{ll}\mathrm{0.172}{r}^{\mathrm{1.28}}{g}^{\mathrm{0.76}}{\mathit{\nu }}^{-\mathrm{0.56}}& \mathrm{if}r\le \mathrm{0.82}\mathrm{mm}\\ \mathrm{0.25}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}& \mathrm{if}r>\mathrm{0.82}\mathrm{mm}\end{array}\right..\end{array}$$

• Step B. During the time T_B, any bubble is able to exchange gas with the surrounding water as a function of u_B and the Reynolds number ${\mathit{Re}}_{\mathrm{B}}=\mathrm{2}{u}_{\mathrm{B}}r/\mathit{\nu }$. In particular, the velocity of gas exchange, j, through a bubble is given by the following expression (valid for Re_B>10):

  $$\begin{array}{}\text{(7)}& j=\sqrt{\left(\mathrm{1}-{\displaystyle \frac{\mathrm{2.89}}{\sqrt{{\mathit{Re}}_{\mathrm{B}}}}}\right){\displaystyle \frac{\mathrm{2}D{u}_{\mathrm{B}}}{\mathit{\pi }r}}}.\end{array}$$

  The velocity of gas exchange is functional to the computation of the gas equilibration time constant of the bubble, T_g, which defines the timescale of the gas transfer process and is given by

  $$\begin{array}{}\text{(8)}& T_{\mathrm{g}}={\displaystyle \frac{r}{\mathrm{3}j\mathit{\beta }}},\end{array}$$

  where β is the Ostwald solubility coefficient, which is the ratio between the volume of absorbed gas and the volume of absorbing liquid for a given condition of pressure and temperature.

• Step C. The gas transfer velocity induced by bubbles, k_B, can be calculated assuming a first-order reaction model which accounts for the time-dependent changes in gas concentration within the drifting bubble:

  $$\begin{array}{}\text{(9)}& k_{\mathrm{B}}={\displaystyle \frac{q_{\mathrm{B}}}{\mathit{\beta }}}\left(\mathrm{1}-e^{-\frac{T_{\mathrm{B}}}{T_{\mathrm{g}}}}\right),\end{array}$$

  with q_B the bubble flux per unit area.

In this application, k_B was calculated by combining the above equations, as detailed in what follows. The bubble rise velocity was calculated assuming a radius r of 0.5 mm using Eq. (6) (leading to u_B=0.133 m s⁻¹). Then, the bubble lifetime T_B was calculated setting α_B to 1.0 in Eq. (5) (T_B≈0.83 s), and the Reynolds number was computed based on its definition (Re_B=133). After that, assuming that in the case of CO₂ $D=\mathrm{1.6}\times {\mathrm{10}}^{-\mathrm{9}}$ m² s⁻¹ and β=0.94, the bubble's transfer velocity, j, and the equilibration time constant, T_g, were estimated via Eqs. (7) and (8) (j≈0.5 mm s⁻¹ and T_g=0.34 s). Finally, k_B was obtained from Eq. (9) using the local flux q_B computed by the numerical code. Given the uncertainty in the value of T_B in this particular experiment, a range of values for the gas transfer velocity induced by bubbles was derived by varying α_B in the interval (0.5, 2.0).

Figure 2 Step-and-pool hydrodynamics. (a) Free surface computed by the FLOW-3D simulation vs. that measured by the hydrometric rod (gray circles) and field of the module of the time-averaged velocity, $\left|\mathit{U}\right|$, along the middle axis of the pool. (b) Time-averaged longitudinal velocity profile, 〈u_x〉(z), in three sections of the pool; the solid blue line with a cyan area denotes the numerical solution, and circles and diamonds denote the 3D and 2D ADV measures, respectively (for the sake of clarity, ADV measures are reported every 10 mm). (c) Streamlines in the first portion of the pool colored according to $\left|\mathit{U}\right|$. (d) Bubbles observed in the pool (blue circles).

3 Results

3.1 Performance of the hydrodynamic model

To assess the reliability of the hydrodynamic model, the depths and velocities measured in the pool were compared to the numerical solution provided by the model. Figure 2a shows the comparison between the observed and modeled water levels along the middle axis of the flume. The water stage ranged from 10.5 in the pool to 12.0 cm downstream of the water jet, whereas the backwater upstream of the plunge had higher water depths. The modeled free surface fit the experimental levels in the pool; over the broad-crested weir; and, with lower accuracy, on the cascade. Therein, however, the measures taken with the rod were sub-optimal owing to observed abrupt variations in the water level. Further, Fig. 2a shows the model's ability to reproduce the module of the velocity vector, $\left|\mathit{U}\right|$. The results highlight a non-uniform flow distribution along the vertical direction within the pool and downstream of the cascade.

Figure 2b reports the profile of the time-averaged longitudinal velocity 〈u_x〉 in three sections located at 10 (A), 20 (B) and 30 cm (C) downstream of the plunging jet. The velocities estimated by the ADV campaign were rather scattered, likely because of the strong fluctuations induced by the turbulent vortexes, and systematically underestimated close to the bottom. Nevertheless, the observed trends were consistent with the numerical solution of the hydrodynamic model. In all cases, the maximum values of 〈u_x〉 were close to the bottom and swiftly decreased downstream of the cascade (from approximately 0.75 m s^−1] (A) to 0.45 m s⁻¹ (C)). Moreover, the 〈u_x〉 decreased in the vertical direction, with a null velocity at 0.45 z₀ for Sect. A and 0.50–0.55 z₀ for Sects. B and C. Negative values are observed at a higher vertical distance from the channel bed, z. Consequently, the maximum backflow was observed near the free surface.

The three examples of velocity profiles pointed out the persistence of a large, counterclockwise eddy structure. The jet plunging in the pool also formed a clockwise vortex in the backwater region, as clearly indicated by the streamlines reported in Fig. 2c. During the experiment, air bubbles entrained by the water jet were visible within a 30 cm wide portion of the counterclockwise eddy and in the backwater region (blue circles in Fig. 2d). The sizes of most of the bubbles highlighted in the picture were in the range r=0.5–1.0 mm, while it was hardly feasible to detect smaller bubbles.

3.2 Gas evasion produced by the step-and-pool configuration

The simulated spatial patterns of total turbulent kinetic energy dissipation rate were highly heterogeneous (Fig. 3a). As expected, the water spilled from the weir dissipated almost all its energy within the region of the two large eddies. However, the highest values of ε were found in the area where the water jet plunges. Therein, ε was higher than 2.0 m² s⁻³ at the free surface, with much smaller values – of the order of 0.1 m² s⁻³ – near the bottom of the flume. The decrease in ε with depth was superlinear (see inset of Fig. 3a). Therefore, overall, the first centimeters of the water column dissipated a significant amount of the jet energy. Moreover, the simulated patterns of ε indicated that, in a small superficial region of the fluid mass with a size of 120 cm² (which is comparable with the footprint of the cascade in the pool), the energy dissipation rate was 1 order and at least 2 orders of magnitude larger than that estimated in the backwater region and downstream of the cascade, respectively.

Figure 3 Scalar fields in the pool. (a) Planar (top) and lateral (bottom) view of the turbulent kinetic energy dissipation rate, ε. (b) Lateral view of the air concentration. (c) Top view of outgas flux, q_B, at three simulation times. Lateral views are along the middle axis of the pool. The right inset in (a) shows the computed ε along z corresponding to the plunging water jet.

The recirculating zones were also able to exchange gas with the overlying atmosphere as they trapped the air dragged by the falling jet (see Fig. 3b). The two large eddies confined the air in some low-velocity regions located in the middle of the water column and swept the bubbles near the bottom and the free surface. According to the numerical model, the maximum bubble concentration was found close to the plunging jet (C≈50 mg L⁻¹), and only a small number of bubbles were predicted at the edges of the region of the observed bubbles. Despite the steady-flow condition achieved (the global air mass in the pool remained almost constant during the simulation), the air entrainment spatial distribution varied significantly over time (see the Video supplement). The magnitude of the bubble flux, q_B, mirrored the patterns of air concentration (Fig. 3c); q_B was highly heterogeneous in the pool, and most air bubbles escaped by bursts in big spots within the first 10 cm downstream of the cascade and the first 5 cm of the backwater region. In contrast with ε, the maximum value (q_B≈5.0 mm s⁻¹) was not found at the cascade, where the air was entrained, but a few centimeters downstream of it.

Figure 4 Gas exchange velocity in the pool. (a) Field of k_T and transverse-averaged exchange velocity ${\stackrel{\mathrm{‾}}{k}}_{\mathrm{T}}$. The solid line is estimated by Eq. (3) with α_T=0.4, while the range of values (colored area) is determined by assuming α_T to be equal to 0.2 and 0.6. (b) Field of k_B and transverse-averaged exchange velocity ${\stackrel{\mathrm{‾}}{k}}_{\mathrm{B}}$. The solid line is estimated by Eq. (9) with α_B=1, while the range of values (colored area) is determined by assuming α_B to be equal to 0.5 and 2.0 (c) The $\mathit{\chi }={\stackrel{\mathrm{‾}}{k}}_{\mathrm{B}}/{\stackrel{\mathrm{‾}}{k}}_{\mathrm{T}}$ ratio within the pool.

The gas evasion induced by the step resulted from the superposition of the effect of turbulence and entrained bubbles. Figure 4a shows the spatial distribution of k_T. We also averaged k_T along y direction (${\stackrel{\mathrm{‾}}{k}}_{\mathrm{T}}$) and reported its longitudinal pattern along the step-and-pool configuration. Owing to the close connection between the gas exchange velocity and ε, ${\stackrel{\mathrm{‾}}{k}}_{\mathrm{T}}$ reaches a maximum of 50 m d⁻¹ near the plunging jet. In this region of enhanced dissipation, local k_T values peaked at 100 m d⁻¹; such a value was almost 1 order of magnitude larger than those observed downstream of the pool, where k_T progressively decreased moving away from the jet. Conversely, in the backwater region k_T was nearly uniform, and ${\stackrel{\mathrm{‾}}{k}}_{\mathrm{T}}$ did not exceed 20 m d⁻¹.

In the region of the domain in which the bubble flow was enhanced, we estimated k_B values which were almost an order of magnitude higher than those computed for the turbulence-driven gas transfer velocity, k_T (Fig. 4b). In particular, k_B locally peaked at 500 m d⁻¹, while the highest values of ${\stackrel{\mathrm{‾}}{k}}_{\mathrm{B}}$ exceeded 100 m d⁻¹ – about twice the maximum value of ${\stackrel{\mathrm{‾}}{k}}_{\mathrm{T}}$.

The patterns of k_T and k_B were not strongly correlated in space and both showed pronounced peaks. Consequently, the relative contribution of bubbles and turbulence to the total gas evasion showed a heterogeneous spatial trend. Figure 4c reports the pattern of the dominance ratio $\mathit{\chi }={\stackrel{\mathrm{‾}}{k}}_{\mathrm{B}}/{\stackrel{\mathrm{‾}}{k}}_{\mathrm{T}}$ along the pool where bubble entrainment was observed. The plot indicates that bubbles dominated the exchange rate (χ>1) in the central portion of the bubble region. In contrast, turbulence led to the outgassing (χ<1) near the cascade (where ε was maximum) and at the edges of the pool (where air concentration was low).

4 Discussion

The important contribution of headwater streams to the global emissions of climatically relevant gases is modulated by the topography of the riverbed, which determines the steep, tortuous and irregular nature of water flow paths in high-energy settings (Duvert et al., 2018; Wallin et al., 2011; Marx et al., 2017). For example, small-scale heterogeneity of channel forms promotes rippled air–water interfaces that enhance the evasion of greenhouse gases. The sequences of step-and-pool configurations typical of mountain streams are plastic examples of small-scale geomorphic elements of rivers that control gas exchange processes at the water–air interface, eventually bearing a relevant impact on the chemical equilibrium of the atmosphere (Botter et al., 2022). Sharp discontinuities of the riverbed of the types observed in steps and small cascades cause a sudden transfer of potential energy to kinetic energy of the water flow. The water jet impacting the pool or the downstream portion of the riverbed then loses a large amount of its energy, roughly corresponding to the height difference between the water levels upstream and downstream of the cascade (i.e., ΔE≈g ΔH, where ΔE is the dissipated energy per unit of fluid mass; see Botter et al., 2022). The dissipation of energy occurs within a few centimeters – a very short spatial scale as compared to the characteristic length of a uniform stream reach with the same energy losses (1–10 m).

In the pool located downstream of a step, the loss ΔE occurs locally through the production and the subsequent dissipation of turbulent kinetic energy. Our experiment highlights that the ε is strongly nonlinear and rapidly decreases downstream of the plunging jet and towards the bottom of the pool. The overall result is that the turbulence depletes indicatively the 15 %–20 % of ΔE near the surface within a length that does not exceed a few centimeters near the cascade.

In addition, the velocity increase observed within the falling jet is responsible for the entrapment and entrainment of the air that overlies the pool. The air bubbles in the flume experiment extended for almost 30 cm downstream of the jet. However, the numerical simulation indicates that most of the air is exchanged only in the first 10 cm near the jet.

In this area, the overall transfer velocity $k=k_{\mathrm{T}}+k_{\mathrm{B}}$ reaches a local maximum of about 500 m d⁻¹, with transverse-averaged values of the mass transfer rate equal to $\stackrel{\mathrm{‾}}{k}=\mathrm{120}$ m d⁻¹. While we recognize that such values of k are quite high – especially in light of the uniformity of the geometry in our flume – we expect that the outgassing induced by steps in mountain streams could be even larger owing to the higher values of ΔH (Natchimuthu et al., 2017; Schneider et al., 2020; Botter et al., 2022) and the enhanced heterogeneity of the bed geometry, which arguably intensify turbulence and air entrainment (Vallé and Pasternack, 2006; Vautier et al., 2020).

Interestingly, in the bubble region of the pool, the potential outgassing, i.e., the surface integral of Eq. (1) for a constant and unit air–water concentration difference, induced by bubble-mediated transport is 2.25 times higher than the corresponding outgassing implied by turbulence alone. Some caution is needed when transposing this numerical result to real-world settings, as the outgassing driven by bubbles may depend not only on the step geometry but also the size and distribution of the bubbles (see Eqs. 6–8). In our simulation, we used a fixed bubble radius, which was set to represent the size of the visible bubbles observed during the experiment. Thus, we have neglected the role of the microbubbles in light of the fact that the mass fraction associated with small bubbles is expected to be negligible (Garrett et al., 2000). Although this analysis is based on several simplifying assumptions, we believe that the orders of magnitude of the underlying processes are properly captured by the model, as space-time heterogeneity in bubble size should not significantly impact the estimated flux of entrained air. Further, we suggest that this result could also apply to different settings, as white waters are observed not only in local steps but also in turbulent reaches where high Froude numbers, macroroughness and emerging stones can promote diffusive entrainment of air. Accordingly, we propose that most of the observed gas evasion in high-energy streams is likely bubble-mediated, in line with previous experimental and modeling studies (Ulseth et al., 2019; Vautier et al., 2020; Klaus et al., 2022). Therefore, we suggest that more effort should be directed to adequately modeling entrained air and gas transfer through bubbles in high-energy streams, where channel bottom morphology causes ripple-free surfaces.

Moreover, the key role of the bubbles in gas evasion from high-energy streams suggested by our results raises two issues about the correct procedure to define the value of k for a given stream reach. First of all, our numerical simulation indicates that the standardization of k should be done not only in terms of Sc but also in terms of β (and D), as suggested in previous studies (Woolf et al., 2007; Hall and Ulseth, 2020; Klaus et al., 2022; Hall and Madinger, 2018). Consequently, some caution is necessary when gas evasion measurements are used under the assumption of k≈k_T if tracer gases with different solubilities are employed. For instance, if in the present analysis low-solubility gases (such as O₂ and He) were considered, the ensuing k_B values would be, respectively, 3.0 and 4.3 times higher than that estimated for CO₂ (${\mathit{\beta }}_{{\mathrm{O}}_{\mathrm{2}}}=\mathrm{0.034}$ and β_He=0.094). Therefore, the effect on the bubble-mediated gas exchange induced by changes in solubility is much higher than that captured by a simple Sc scaling.

Second, in spite of the observed empirical correlation between k and $\stackrel{\mathrm{‾}}{\mathit{\epsilon }}$ in high-energy streams (Ulseth et al., 2019), our results indicate that the definition of a robust upscaling law between the gas transfer rate and the turbulent kinetic energy dissipation rate needs to include the spatial frequency of steps and the relative contribution of bubble-mediated and turbulence-induced gas evasion, as long as the scaling of k_B could be modulated by local hydrodynamic features such as the flow velocity, the water depth and the channel morphology.

A final remark concerns the local nature of the outgassing process in our step-and-pool configuration and the related implications for future gas transfer studies. Here, 75 % of the bubble-mediated potential outgassing occurs through spots with high q_B, the area of which does not exceed 2.5 dm² overall – less than 20 % of the region where air entrapment is visible. This heterogeneity of gas evasion in the pool is also mirrored by the spatial patterns of k_T. The first 5 cm downstream of the cascade caused almost 25 % of the potential outgassing induced by turbulence in the whole domain. Thus, both the turbulence-induced and bubble-mediated processes mostly take place close to the plunging jet and operate at a scale that is much smaller than the length used to measure k during in-field experiments, which ordinarily adopt probes or chambers with an inter-spacing of at least some meters (Vingiani et al., 2021; Vautier et al., 2020; Botter et al., 2022). This difference in scale might explain the heterogeneity of k that can be empirically estimated from tracer gas concentration drops observed along reaches or segments that contain steps or cascades (Wallin et al., 2011; Natchimuthu et al., 2017; Leibowitz et al., 2017). In these instances, k is implicitly averaged over the entire fluid region included within the reference sampling points and cannot capture the local peaks of k, as the step contribution to the outgassing is eventually “smoothed” over a larger area. An important consequence of the local and heterogeneous nature of the outgassing process in steps and pools is that, in cases where the spatial patterns of local k are unknown, the measure of the transfer velocity is dependent on the position of the monitoring points and the geometry of the steps. Thus, the average value of k across a step is not only unpredictable but also inevitably scale-dependent (i.e., the larger the water volume involved in the averaging, the lower the effective mass transfer velocity k; Botter et al., 2021, 2022). As a consequence, aggregating k across heterogeneous reaches might be inherently biased, and extrapolating observed reach-wise mass transfer rates could be appropriate only in a limited subset of reaches that share the same degree of heterogeneity of the reach where the experiments were performed. Based on the above arguments, we propose that the use of the mass transfer rate, k, should be used with caution when the heterogeneity of the flow field controls the fraction of mass evaded into the atmosphere, as in our step-and-pool configuration. In these cases, the mass transfer rate could be highly site-specific, and its extrapolation to different settings could be feasible only following a detailed analysis of the evasion produced by the constitutive geomorphic elements of the focus reach (including steps and cascades) (Botter et al., 2022).

5 Conclusions

In this study, we have numerically simulated the hydrodynamics of a step-and-pool configuration, which were also reproduced in a laboratory-designed flume experiment. Computational values of key hydrodynamic parameters that control the gas transfer between water and the atmosphere were then used to estimate the spatial patterns of gas transfer velocity and the contributions to the gas evasion due to turbulence and air entrainment. Our simulation shows that the cascade downstream of the step dissipates energy and entraps air in a small, irregular region near the chute (about 3 dm²), in which the mass transfer velocity k is locally very high (500 m d⁻¹). This result indicates that gas evasion in step-and-pool configurations may be a very local process, taking place within a few square decimeters of the plunging jet. The numerical simulation also suggests that bubble-mediated gas exchange dominates the turbulence-induced outgassing, reinforcing the idea that gas evasion in mountain streams could be mainly driven by bubbles and white waters. The observed heterogeneity of k – the pattern of which is linked to the specific geomorphic features of the step under investigation – raises concerns about the ability of traditional metrics such as the mass transfer rate to quantify gas evasion from step-and-pool configurations and contrast the outgassing potential of different morphologically heterogeneous streams.