Canopy Randomness, Scale, and Stem Size Effects on the Interfacial Transfer Process in Vegetated Flows

Aquatic vegetation plays an important role in natural water environments by interacting with the flow and generating turbulence that affects the air‐water and sediment‐water interfacial transfer. Regular and staggered arrays are often set as simplified layouts for vegetation canopy to study both mean flow and turbulence statistics in vegetated flows, which creates uniform spacing between vegetation elements, resulting in preferential flow paths within the array. Such preferential paths can produce local high velocity and strong turbulence, which do not necessarily happen in natural environments where vegetation is randomly distributed. How the randomness of the canopy affects interfacial processes by altering spatial turbulence distribution, which can potentially lead to different turbulence feedback on the interfacial transfer process, remains an open question. This study conducted a series of laboratory experiments in a race‐track flume using rigid cylinders as plant surrogates. Mean and turbulent flow statistics were characterized by horizontal‐ and vertical‐sliced PIV. Based on the measured flow characteristics under different stem diameters and array configurations, we propose a method to quantify the randomness of the vegetation array and update a sediment‐water‐air interfacial gas transfer model with the randomness parameter to improve its accuracy. The updated model agrees well with the dissolved oxygen experimental data from our study and data from existing literature at various scales. The study provides critical insight into water quality management in vegetated channels with improved dissolved oxygen predictions considering vegetation layout as part of the interfacial transfer model.

Regular and staggered arrays are often set as simplified layouts for vegetation arrays to study the mean flow and turbulence statistics in vegetated flows (Dupuis et al., 2016;Etminan et al., 2018;Li et al., 2020;D. Liu et al., 2008;Lopez & Garcia, 2001;Stoesser et al., 2010;Yang et al., 2015), which creates uniform spacing between vegetation elements, resulting in preferential flow paths within the array.Such preferential paths could produce local high velocity and strong turbulence (Chang & Constantinescu, 2015;Koken & Constantinescu, 2021), which does not generally happen in natural environments where vegetation is randomly distributed (Lara-Romero et al., 2016).Hence, some studies have opted to set up random arrays as the canopy layout (Nepf, 1999;Tanino & Nepf, 2008a;Tinoco & Coco, 2016;Tinoco & Cowen, 2013).Etminan et al. (2018) showed that the "sheltering effect (Nepf, 1999)" between two in-line neighboring stems could lead to a large difference of drag predictions on individual stems.Such an effect is more prominent when the array is in random distribution compared to the staggered configuration, due to the lower incident velocity faced by the downstream stems compared to the upstream ones.The "blockage effect," on the other hand, addresses the local increase in the velocity around the vegetation due to the presence of the stem elements reducing the cross-sectional area locally, which is expected to be less prominent in a random canopy because of less preferential paths, created by the uniform spaces between vegetation elements, existing inside the canopy compared to a regularly distributed layout.Shan et al. (2019) further investigated the differences in drag distribution between random and aligned vegetation.By comparing the experiment results in uniform and clumped, heterogeneous arrays, Park and Hwang (2019) indicated that under the same canopy density, the longitudinal dispersion coefficient becomes higher as the vegetation pattern becomes more random, which could also affect sediment suspension and bedload transport rate (Shan et al., 2020).However, only a few studies have focused on the effects of spatial heterogeneity and the resulting hydrodynamics.Ricardo et al. (2016) conducted laboratory experiments to investigate the effects of spatial nonuniformity of vegetation arrays by varying the longitudinal vegetation density.The results revealed that the time-averaged mean flow field is highly influenced by the locally varying vegetation density due to the upstream stem distribution as a form of spatial memory (see the conceptual sketch in Figure 1b).Tanino and Nepf (2008b) developed a model to predict the lateral dispersion in random vegetation arrays, considering the spatially heterogeneous velocity field, which arises from the random distribution of the stem elements.However, there is no well-recognized approach to quantify the randomness level of a plant canopy based on the measured flow characteristics, which could lead to different representative turbulent quantities.Tseng andTinoco (2020, 2022) addressed that the vegetation-generated turbulence plays an important role in the AWI and SWI transfer process (see the conceptual sketch in Figure 1a), but how randomness of the canopy and the size of vegetation stems affect the transfer process by altering the turbulence flow pattern remains an open question.
Two major contributions of this work are: (a) to provide a method to quantify the randomness of the vegetation array based on the measured flow characteristics, and (b) to investigate the randomness and stem-size effects on the interfacial gas transfer processes across AWI and SWI, to improve predictions by models proposed by Tseng andTinoco (2020, 2022).This study, for the first time, accounts for vegetation distribution patterns in the interfacial gas transfer models to improve predictions for more universal applications to the aquatic ecosystem environments.The article is structured as follows.In Section 2, we describe the methodology, experimental setup, and measurements of flow hydrodynamics and DO.In Section 3, we present the horizontal and vertical-sliced PIV results of representative mean and turbulent flow quantities under different stem diameters and array configurations according to the submergence ratios, and provide a quantification method to discuss the randomness level of the array.We then show the DO measurement results for the AWI gas transfer rate and the SWI diffusion coefficient.To conclude, we introduce the newly proposed canopy randomness indicator into the prediction models for AWI and SWI interfacial gas transfer and compare the predicted results to experimental data.

Experimental Setup
A series of laboratory experiments were conducted in a race-track flume (RTF) (Figures 2a and 2d) at the Ecohydraulics and Ecomorphodynamics Laboratory at the University of Illinois at Urbana-Champaign.Flow is driven by a paddle wheel, generating uni-directional, recirculating flows, which is controlled by an inverter whose frequency, f, was set in a range between 2 and 6 Hz.The rotation speed of the paddle wheel, ω, can be estimated based on a linear relation form: ω(rpm) = 0.79 f(Hz).The chosen frequencies resulted in mean bulk velocities, U = 0.05-0.18m/s, yielding flows with Reynolds number, Re H = UH/ν, where H is the water depth and ν is kinematic viscosity, ranging from 3,000 to 22,000, covering a range of sediment motion from no suspension to suspension.
The recirculation feature of the flume provides an infinite length for flow to develop (Piedra-Cueva et al., 1997), while the paddle wheel is designed to produce a uniform spanwise velocity profile with minimal disturbance to ensure full development of the flow and turbulence features while continuously recirculating.The 0.75 m channel width minimizes the effects from the side walls.Bends were equally separated into four waterways by three baffles to minimize overturning currents and consequent energy losses, and the 9.57 m long straight test section could further break down secondary flow structures generated in the bending regions.
Well-sorted crushed walnut shells, with a median grain size, D s ≈ 1 mm, were used to construct a 10 cm thick sediment bed.Compared to the typical sand whose specific gravity ≈2.65, crushed walnut shell is relatively light with a low specific gravity ≈1.48 (measured according to ASTM ( 2014)), making it an ideal sediment substrate for scaling sediment dynamics in laboratory-scale experiments (Heller, 2011;Kleinhans et al., 2014).The settling velocity, w s , is approximately 1 cm/s estimated by Rubey (1933).
To compare with previous studies on the interfacial gas transfer process in vegetated flows (Tseng & Tinoco, 2020, 2022), we followed their setup for the vegetation surrogate.Arrays of rigid acrylic cylinders with height, h = 10 cm, were used to simulate different sizes of high-stiffness aquatic vegetation (e.g., Sagittaria sagittifolia and Sparganium erectum) using two diameters, d = 0.635 and 2.54 cm.Two layouts, staggered and random (Figure 2c), of the uniformly distributed array under the same array density, were designed to investigate the effects of different vegetation distribution patterns.The solid volume fraction, ϕ, was set as 0.5 for both cases with different stem sizes, while roughness density, ah, reduces from 1.0 to 0.25 when d increases from 0.635 to 2.54 cm.For the staggered array, we set the streamwise distance between each row as Δx = 2d, and the spanwise distance between each column as Δy = 4d under a staggered configuration.For the random array, we followed the methodology described in Tinoco and Coco (2016), using a MATLAB (MathWorks, Inc) routine to randomly pick the populated locations for each element under the same domain size as the staggered array with a minimum Water Resources Research 10.1029/2023WR035220 allowed distance between adjacent elements (s min = d).Once all the locations have been selected, we plotted the changes in density across both streamwise (x) and spanwise (y) directions to ensure a uniform distribution across the whole area.If not, then we repeat the same procedures until a satisfactory uniformity is reached.The total length of the vegetation canopy is approximately 5 m, starting at x = 2.57 m, and the center of the measurement section was set at x = 6.38 m to allow the flow to develop and avoid influence from the upstream bend (Figures 2a  and 2b).Two submergence ratios were used: h/H = 0.5 (half-submerged) and h/H = 1 (emergent).
There were 24 runs with 3 different flow velocities conducted in the RTF, each under 2 submergence ratios (h/ H = 0.5 and h/H = 1), 2 stem diameters (d = 0.635 and 2.54 cm), and 2 array configurations (staggered and random).Detailed physical parameters are shown in Table 1 below.

Hydrodynamic Measurements
Instantaneous velocities were measured by 2D Particle Image Velocimetry (PIV), focusing on a small gap region (0.75S x for sparse and 3S x for dense cases, where S x is the mean distance between vegetation stems in x-direction in staggered arrays) as shown in the red dashed area in Figure 2b.The gap was left to avoid laser reflection issues by the dense acrylic cylinders within the array.As discussed in Tseng and Tinoco (2022), the measurement uncertainties inside the gap compared to the measurement inside the canopy depend on the choice of the analysis window, which can be reduced by analyzing the area within the gap in dimensions similar to the spacing within the array.Hence, to ensure representative measurements and consistency throughout all the data sets, the lengths of the window for analysis were set as 4d for all cases (blue dashed area in Figure 2b), the same as S x , which is small enough to avoid significant gap issues.
A 5 Watt (W) Continuous-Wave 532 nanometer (nm) laser with a cylindrical lens was used to generate vertical and horizontal planar light sheets for PIV, with a thickness of <1 mm along the center row, y = 0, of the array for the vertical slice and at the height z = 0.5 hr, for the horizontal slice.Neutrally buoyant particles (silver-coated hollow glass spheres) of 10 micrometers (μm) in diameter were used as PIV seeding.A 5-Megapixel CCD Camera, JAI GO-5000M-USB3, with a Navitar 25 mm focal length lens was used to capture 8-bit grayscale images for 1 min in each run, and the recording frame rate was set between 40 and 60 Hz depending on the flow speed.A convergence analysis presented in Text S1 in Supporting Information S1 shows reasonable convergence over this 1 min recording setup.A MATLAB-based open-source tool, PIVlab (Thielicke & Stamhuis, 2014), was used to process raw images with their corresponding spatial calibration factor, CF (cm/pixel).CF goes from 0.009 to 0.017 for vertical PIV and 0.025 for horizontal PIV.Three consecutive passes with a 50% size reduction between passes and 50% overlapped interrogation areas were used to obtain higher resolution results in crosscorrelation calculations, and the subwindow size of the first pass followed the one-quarter rule (Keane & Adrian, 1990).The post-processed instantaneous velocity components in streamwise x-, lateral y-, and vertical zdirections are denoted as u, v, and w, respectively, which can be used to calculate the velocity fluctuations (u′, v′, w′) using Reynolds decomposition as: where the operator, "〈 〉," represents time averaging, following the notation used in Tseng andTinoco (2020, 2022).

Air-Water Surface Gas Transfer Rate
In river flow systems with a permeable sediment bed, gas can transfer across surface and bottom interfaces.The depth-integrated advection-diffusion equation of C DO , the spatial-averaged DO concentration away from the bed (see Figure 1a), under full developments in x-and y-directions can be used to describe the aeration process with a DO source across AWI, J AWI , and a DO sink across SWI, J SWI : In natural rivers, flows are generally shallow and fast.Under this condition, J SWI is often much smaller than J AWI (e.g., Hartman & Hammond, 1984), Equation 2is then mainly controlled by the aeration flux from the air: According to the methodology proposed by the American Society of Civil Engineers ( (Stenstrom, 2007)), Cobalt Chloride Hexahydrate (CoCl 2 ⋅ 6H 2 O) was put into water as a catalyst to mix with Sodium Sulfite (Na 2 SO 3 ) as the oxygen depletion agent to deplete all the DO in water.Once DO concentration reached the minimum level close to zero, fresh oxygen gas would be re-aerated from the surface.The whole re-aeration process was monitored by two DO sensors (PASPORT Optical Dissolved Oxygen Sensor, PASCO) at sampling intervals from 30 to 60 s, located near the surface and near the bed to monitor bulk and near-bed DO concentrations.Based on the Diffusive Boundary Layer theory (Danckwerts, 1951) with Fick's diffusion law: where C w is the local instantaneous DO concentration in water (see Figure 1a), the recovery curve of DO concentration under different scenarios can be used to fit the corresponding air-water surface gas transfer rates, k L , by solving the above first order ordinary differential equation with an assumption that C DO t=0 = 0: where H is the averaged water depth and C sat is the saturated DO concentration.The above re-aeration equation was widely used for experimental designs to measure the gas transfer rate, k L , from the re-aeration process (Leman et al., 2018;Stenstrom, 2007;Tseng & Tinoco, 2020, 2022).

Diffusion Coefficient Across Sediment-Water Interface
Gas transfer across SWI depends on both the bio-chemical consumption inside the sediment layer and physical diffusion across the interface.In a small-scale flume experiment with a shallow and fast flow, the physical diffusion process dominates the SWI gas transfer as we can see that ΔC DO , the difference between C DO and the near-bed spatial-averaged DO concentration, C DO,b (see Figure 1a), approaches zero at the end of the re-aeration process (see Figure 6b in Tseng and Tinoco (2022)).However, in a large-scale flume experiment with a relatively slower flow system according to its domain scale, Figure 3 shows that ΔC DO approaches a constant higher than zero at the end of the re-aeration process, indicating that some diffusible reductants, such as Fe ++ and Mn ++ (Bouldin, 1968), were growing during re-aeration in the sediment layer.The reductants kept consuming DO in the bed, which caused a ΔC DO , driving physical diffusion across SWI.According to Figure 3b, ΔC DO increases from a lower level at the beginning of re-aeration and reaches a higher equilibrium level at the end.Following the first-order reaction model proposed by Bouldin (1968), the bio-chemical consumption rate, J bio chem , is proportional to C DO,b : where k sed is the first-order rate constant, which can be estimated according to the final equilibrium state: t eq is the time when the system reaches equilibrium.

Water Resources Research
10.1029/2023WR035220 The diffusive flux at the SWI consists of three components: molecular diffusion, dispersion, and turbulent mixing (O'Connor & Harvey, 2008;Voermans et al., 2017Voermans et al., , 2018)): where θ is the sediment porosity, z = 0 denotes the vertical location at the SWI, and D eff is the effective diffusion coefficient, which is the sum of three diffusion coefficients: the molecular diffusion coefficient, D, the dispersion coefficient, D d , and the turbulent diffusion coefficient, D t .The value of J SWI depends on both the effective diffusion coefficient, D eff , and the near-bed DO concentration gradient, ∂C DO ∂z | z=0 .The theoretical DO profile shown in Figure 1a suggests that DO concentration quickly reaches a constant above the diffusive layer all the way up to the surface, which approaches the value of C DO measured near the surface.Then the near-bed DO concentration gradient can be approximated as by assuming that the DO concentration profile is linear within the diffusive layer.Replacing Equation 9 back into Equations 8 and 7, we can thus calculate k sed in each case: δ c is the thickness of the diffusive layer, which can be estimated by using the eddy viscosity distribution model at the top diffusive boundary layer (Dade, 1993) according to Reynolds analogy: where u* is the bed shear velocity; β is an empirical coefficient that can change with sediment size and permeability.However, previous experiments suggested that β = 19.4 ± 5.5 (Steinberger & Hondzo, 1999) with fine sediment and low permeable bed, which holds the same order of magnitude in a certain range of the permeability Reynolds number, Re k (= u * ̅̅̅ ̅ K √ /ν, where K is the sediment permeability).In vegetated flows, following Tseng and Tinoco (2022), instead of using u*, δ c must be corrected by the total near-bed TKE, k t , to consider the vegetation-generated turbulence: By further integrating with time through the whole re-aeration process, the integral of the SWI flux, J SWI , is equal to the integral of the bio-chemical consumption flux, J bio che , plus the total saturated DO stored inside the sediment bed per unit area, Q DO,s : where t = 0 is the starting time of re-aeration, T is the re-aeration time period.Finally, the effective diffusion coefficient of SWI flux, D eff , can be obtained by solving the above integral equation.

Hydrodynamic Features
In uni-directional canopy flows, the magnitude of streamwise velocity (u) is much generally larger than the spanwise and vertical velocities (v and w).However, when doing streamwise average analysis to investigate the

Water Resources Research
10.1029/2023WR035220 canopy-flow features, the averaging process will effectively smooth out the streamwise variations (∂/∂x) and present just the mean value for the property.At this condition, only the spanwise and vertical variations (∂/∂y and ∂/∂z) are important.As emergent and submerged vegetation show completely different mean and turbulent flow features from the streamwise average analysis, two dominant components of turbulence are vertical shear turbulence and horizontal wake turbulence (Nepf, 1999;Nepf & Vivoni, 2000) that was measured by the verticalsliced PIV and horizontal-sliced PIV, respectively.TKE production by the two components of turbulence has been shown as the two representative indicators for air-water surface gas transfer prediction, depending on the submergence conditions (Tseng & Tinoco, 2020, 2022).
For submerged cases (h/H = 0.5), Figure 4 shows the normalized mean flow field, mean flow velocity profile, TKE profile, and TKE production profile measured from the vertical-sliced PIV at the center row of the vegetation array.Since the vertical-sliced PIV can only measure the streamwise (u) and vertical (w) velocity, TKE was approximated as: This approach was justified based on previous measurements by Tanino and Nepf (2007), indicating that u′ 2 ≈ v′ 2 in the canopy flow, which was generally used to approximate the estimation of TKE in 2-D flow measurements (Tanino & Nepf, 2008b;Tseng & Tinoco, 2020;Yang et al., 2016).The vertical-shear TKE production, P uw , was then calculated based on the vertical Reynolds stress, 〈u′w′〉, times the time-averaged vertical velocity gradient: When averaging P uw in streamwise direction across the canopy, we ignored the additional dispersive production, P uw,dis = 〈u′′w′′〉 ∂〈u〉 ∂z , formed from the spatial correlations in the time-averaged velocity field (Lopez & Garcia, 2001) by approximating the streamwise-mean production as the total averaged production, which has been TSENG AND TINOCO proved valid when canopy density is not extremely sparse (ah > 0.1) (King et al., 2012;Poggi et al., 2004;Tseng & Tinoco, 2020).
As shown in Figure 4, the discontinuous drag by the submerged canopy creates a sharp velocity gradient that induces a strong shear flow on top of the vegetation array.The shear generates canopy-scale vortices that rotate relatively to the y-axis, which results in a peak value of TKE and P uw .As discussed by Tseng and Tinoco (2020), these canopy-scale vortices dominate the vertical mixing across the upper-canopy and the within-canopy layers, and can further enhance the surface interfacial gas transfer as predicted by using the peak value of P uw as the indicator.When comparing the results between different layouts (random and staggered) under two stem sizes (thick and thin), profiles are similar between the two canopy layouts when d = 0.635 cm (Figures 4a and 4b).The similarity in these profiles is explained by noticing that when d = 0.635 cm, the array is considered dense (ah = 1.0) (Coceal & Belcher, 2004), and the difference of "sheltering effect" and "blockage effect" (Etminan et al., 2018) is not evident, resulting in similar flow variations along the spanwise direction.In this case, verticalsliced PIV measured along the centerline for both layouts shows equally strong shear from the upstream vegetation observed in the measuring window.
However, when d = 2.54 cm, the array (ah = 0.25) is considered sparse (Coceal & Belcher, 2004), and the spanwise variations in the flow are more different from each other.Profiles under random array show a milder velocity gradient on top of the canopy and weaker peak values of TKE and P uw profiles in the random canopy (Figure 4c), compared to the profiles in the staggered case that show stronger peaks since the PIV slice followed exactly along the center row of the array.The difference caused by the upstream "sheltering effect" and "blockage effect" is now prominent between the two layouts, which indicates that the spanwise variation is important and affects how representative the peak values of TKE and P uw are when using them as key predictors in interfacial gas transfer models.This effect is seemingly more important for relatively sparser arrays and larger individual elements, which will be discussed in more detail in Section 3.3.
In the emergent case, the presence of vegetation stems exerts uniform drag on the flow, making a nearly constant velocity profile above a certain height from the bed (Figure 5).The presence of vegetation redistributes the pressure field and separates the bottom boundary shear layer, lifting it up to form the prevailing coherent vortex structures (e.g., horseshoe vortex) (D. Liu et al., 2008).These temporal-varying vortex structures can entrain high momentum flow from the surroundings of the plants into the wakes as shown by the characteristic velocity bulge near the bed from the velocity profile, which was also presented in previous studies (Dupuis et al., 2016;Stoesser et al., 2010;Tseng & Tinoco, 2020, 2021).As discussed in the submerged cases, due to more prominent differences in the upstream "sheltering effect" and "blockage effect" between two canopy layouts, similar velocity fields under random array and staggered array are presented in Figures 5a and 5b when d = 0.635 cm.However, when d = 2.54 cm, the thicker stems create a wider low-pressure region downstream of each stem, which yields much lower velocity at x/d = 0 right behind a stem and allows for more rapid flow recovery downstream from the vegetation in the staggered array (Figure 5d), which cannot be seen in the random array because of a more prominent "sheltering effect" and less "blockage effect" (Figure 5c).
The centerline vertical-plane PIV measurement shows different features of emergent and submerged vegetation canopy.However, it is important to further understand the horizontal variations of the flow so as to investigate the effect of different stem distributions (e.g., staggered and random).For example, flows along the preferential paths show totally different features from flows along the centerline of the staggered vegetation canopy, creating different horizontal flow variations under different canopy layouts.More detailed discussions are presented in Text S2 in Supporting Information S1.To characterize the spanwise variations in mean and turbulent flow quantities due to the randomness of the array, we conducted PIV on a horizontal slice located at the mid-height of the vegetation (z = 0.5 hr).In the horizontal-sliced PIV, TKE was defined based on u′ and v′, following the same approximation approach used for vertical-sliced PIV in the submerged case (Equation 14): and the horizontal-shear TKE production, P uv , was calculated based on the horizontal Reynolds stress, 〈u′v′〉, times the time-averaged spanwise velocity gradient: When spatial-averaging, we again ignored the dispersive flux since the canopy density is higher than 0.1 in all experiment cases according to Equation 16 where " -" represents spatial average in spanwise direction.Figure 6 shows the normalized mean spanwise flow field, mean spanwise flow velocity profile, TKE profile, and TKE production profile measured from the horizontal-sliced PIV in emergent vegetation canopy.When the canopy is emergent, vegetation covers the entire water column, exerting continuous drag to the flow without the sharp vertical velocity gradient present in the submerged canopy.Under this condition, horizontal stem-scale wake turbulence is dominant (Nepf, 1999;Tanino & Nepf, 2008b).Since stem-generated stem-scale turbulence is in general uniformly distributed along the z-axis as the mean flow velocity profile (Tseng & Tinoco, 2020), mid-height horizontal-sliced PIV data provides representative mean and turbulent flow quantities for the array.
Under staggered configuration, the spanwise velocity profile shows a regular, periodic pattern, where the crosssectional flow variations can be clearly distinguished.However, for the random array case, spanwise velocity profile becomes irregular and non-periodic.Peak values of TKE and P uv are found at the inflection points of the velocity profile, due to horizontal shear gradients that depend on the distance between individual elements.Spanwise TKE and P uv patterns are closely correlated to the velocity profiles, as we can also see regular, periodic profiles in staggered cases and irregular, non-periodic profiles in random cases.This feature is important and being followed for quantifying the randomness of the array, which will be discussed in detail in the next section (Section 3.2).

Randomness Effects
The horizontal-sliced PIV data (Figure 6) clearly presents the spanwise variations of the flow.We applied Fast Fourier Transform (FFT) to convert the discrete spanwise-varying flow data into the corresponding spectrum domain: where N is the total number of spatial elements in the velocity data, A k is the amplitude corresponding to the kth wavenumber in the spectrum, F * n represents the nth flow data point along the y-direction.Following Equation 20, FFT decomposes the spanwise flow pattern into its wavenumber components, which provides information about the periodicity characteristic and pattern regularity for understanding the randomness level.Since TKE and TKE production are in dimensions of 〈u〉 2 and 〈u〉 3 , respectively, we focused on the FFT conversion of the spanwisevarying F * = 〈u〉 2 and F * = 〈u〉 3 as shown in Figure 7, which closely correlates to the regularity of the spanwise mean flow velocity pattern.The spectra is normalized by the total area under the curve, which is the total amplitude contributions among all the wave numbers.Comparing the two array layouts, spectra tend to focus more on higher peak values with a narrower bandwidth in the staggered cases, while peaks are lower and with a wider bandwidth of the spectrum curve in the random cases, especially when d = 2.54 cm.
According to Ranjan et al. (2022), the measuring gap regions affected by the upstream vegetation stems have a short extent.Turbulence strength downstream of the array reaches the maximum value within one unit of the diameter length, 1 d, and then transitions from an exponential to a power-law decay, quickly dropping back to a minimum value within four units of the diameter length, 4 d, suggesting that the closest two rows of the canopy upstream from the measuring gap would have the highest influence to the flow pattern.Hence, we propose a newly defined randomness indicator, ξ, based on the corresponding FFT spectra of 〈u〉 2 or 〈u〉 3 in Figure 7 to quantify the randomness level of the canopy: where A 1st and A 2nd are the first and second highest amplitudes in the spectrum, and the value of ξ is in a range between 0 and 1.When the canopy is perfectly staggered, due to consideration of the very first two upstream rows, the ideal flow pattern would be constructed perfectly by the combination of two sine waves, which leads to ξ = 1 since the amplitude contributions in the spectrum are only from the two dominant wavenumbers.However, since there are always measurement deviations in the experiment depending on the resolution of the measurement, ξ will never be perfectly one.Instead, we can obtain a value that is close to one as ξ ≈ 1 in staggered cases (Figures 7b and 7d).When vegetation distribution becomes more random, the two upstream rows would result in a relatively irregular pattern (Figures 7a and 7c) and result in a lower value of ξ depending on how random the array is since the amplitude contribution in the spectrum is now affected by several wavenumbers in a wide range.In the next section, we implement this ξ value into previously proposed interfacial gas transfer models for AWI and SWI (Tseng & Tinoco, 2020, 2022) to account for the array configuration effects.(Tseng & Tinoco, 2020) proposed a modified SR model to predict air-water surface gas transfer rate, k L , in vegetated flows by using the characteristic shear velocity, u * c , and the bulk-scale characteristic TKE production, P c , to normalize the commonly used SR model for open-channel flows:

Interfacial Gas Transfer: Modified Surface Renewal Model
L + is a non-dimensional length scale consisting of a characteristic Reynolds number, a characteristic shear velocity, and a characteristic velocity scale: where U c is the characteristic mean flow velocity.
For the emergent case, the presence of vertically mounted vegetation stems throughout the entire water column induces horizontal shear, generating stem-scale wake turbulence which dominates the transfer process.Therefore, instead of using PIV data from the vertical slice (Tseng & Tinoco, 2020), horizontal-shear mean TKE production, , and horizontal-mean flow velocity, U m , measured by the horizontal-sliced PIV, were used as P c , u * c , and U c in the prediction of k L as: α is an empirical coefficient that is determined by the experimental data according to the above equation.In this study, we used MATLAB Curve Fitting Tool (cftool) with linear least squares fitting to obtain the best-fitted α.L eme + is the non-dimensional length scale for the emergent vegetation: As shown in Figure 8a, data from both staggered and random cases already collapse and follow the previous data (Tseng & Tinoco, 2022) and the model predictions without doing any randomness factor correction, yielding the best-fitted α = 0.12 when Re d is larger than the critical Re dc = 300.When considering the bulk-mean turbulence quantities as predictors in the SR model, the average turbulence magnitude doesn't change too much under the same array density and flow velocity (see TKE and P uv profiles in Figure 6), as proven by C. Liu et al. (2021) and Zhao and Nepf (2021), addressing that effects of stem size and canopy layout can be negligible on turbulence strength and sediment bedload transport.
In the submerged cases, vertical-shear-induced canopy-scale turbulence on top of the vegetation array dominates the transfer process.According to Tseng and Tinoco (2020), the maximum streamwise-averaged vertical-shear TKE production, P uw,max , vertical-shear velocity, , and vertical-sliced upper-canopymean velocity, U up , measured from the vertical-sliced PIV can be used as the key predictor for the modified SR model: Water Resources Research 10.1029/2023WR035220 TSENG AND TINOCO L sub + is the non-dimensional length scale for the submerged vegetation where Re up = U up L up /ν is the upper-canopy mean flow Reynolds number and L up is the upper-canopy length defined as the canopy to surface distance.However, since P uw,max is measured based on the peak value on top of the vegetation along the center row of the canopy, we need a correction parameter to determine how representative the peak value is for the whole array, which would expect to be highly sensitive to the canopy layout.
Assuming that the spanwise variation pattern of P uw,max generally follows the spanwise variation pattern of 〈u〉 3 , and P uw has dimensions of h〈u〉 3 , we determined a randomness correction parameter, η, which is the ratio of the previously defined array randomness indicator measured under the random layout, ξ rand , to the value measured under the staggered layout, ξ stag , for 〈u〉 3 spectrum with the same stem size and array density: We then introduce η into Equation 26 to consider the representation of P uw,max in different layouts Figure 8b shows that the experiment data under different layouts and stem sizes collapse well into the model prediction line when the upper-canopy Reynolds number is lower than the critical value, Re up < Re upc = 25, 000, yielding α = 0.065.Previous data from (Tseng & Tinoco, 2020, 2022) gives α = 0.71 when Re up > Re upc .
However, the performance of the FFT method is highly affected by the data resolution.Hence, when using the η correction, we need to make sure the measurements in random and staggered arrays have a similar spatial resolution so that the new updated model (Equation 29) can provide improved predictions based on the comparison between the two cases.

Modified SWI Gas Transfer Model
According to the Re TKE -dependence model proposed by Tseng and Tinoco (2022) to predict the effective interfacial gas transfer coefficient, D eff , in vegetated flows, the model considers total near-bed turbulence generated from the bed shear stress, k tb , the vegetation stem, k tv , and the instantaneous non-linear coherent structures, k tcs : The first two terms, k tb and k tv , can be approximated by the empirical relations proposed by Stapleton and Huntley (1995) and Tanino and Nepf (2008b), respectively, which have been widely applied to estimate these two components in vegetated flows (Tseng & Tinoco, 2021;Yang & Nepf, 2018, 2019;Zhao & Nepf, 2021) when vegetation canopy is neither too sparse (ϕ > 0.01) nor too dense (ϕ < 0.3), using: where u* is the bed shear velocity, and C ub is the integral coefficient that converts mean flow velocity to bed shear velocity based on log-law (Wilcock, 1996) κ = 0.41 is the von Karman constant and k s is the bottom roughness, which is often set as the median grain size of the sediment particles, D s .
where C D is the drag coefficient that is approximately 0.8 for ad ≈ 0.06 in this study (Nepf, 1999).The last instantaneous non-linear term, k tcs , can be estimated by a scaling method according to the square root of the combined time-averaged k tb and k tc with a non-linear empirical coefficient, c, that depends on the vegetation shape, canopy array density, and the stem Reynolds number, Re d (Tseng & Tinoco, 2021): For cylinder arrays that have enough space for full development of coherent structures, c is suggested in a range between 4 and 5 (Cheng et al., 2018;Tseng & Tinoco, 2021).Then Re TKE can be calculated based on k t As discussed in Section 3.3, the averaged bulk mean k tv would not be significantly affected by the stem size and the array layout.However, in the randomly distributed canopy, due to the sheltering effect between two adjacent vegetation rows, downstream stems would face a lower incident velocity than the upstream ones, which could lead to a large difference in the time-averaged estimation of the resulting turbulence by the instantaneous fluctuating coherent structures according to Equation 34.Hence, we corrected the empirical non-linear coefficient, c, Equation 34 for k tcs is then updated as: Figure 9a shows the modified Re TKE -dependence model prediction by Tseng and Tinoco (2022)  showing an overestimation of Re TKE when the canopy is dense.This is because, in the dense canopy (ah = 1.0), the limited space between vegetation elements restricts the development of such temporal-varying vortex structures, reducing the strength of these coherent systems closely surrounded by the viscous layer around the stems' surface.Under this condition, the originally downstream-stretching vortices would be broken down into smaller eddies (Pope, 2000) with a weaker component from the k tcs term, suggesting that the coefficient c would need to be updated with a smaller value.Replacing c with a new best-fitted value c = 0.9 for the dense canopy cases (d = 0.635 cm, ah = 1.0), all data from channels with vegetation under different stem sizes and array densities (color markers) now matches better with data from bare-bed channels (gray markers) and follows the updated model predictions as shown in Figure 9b.The total normalized variation was calculated as: which reduces from 4.23 to 1.16 after the coherent vortex correction as shown by the insets in Figures 9a and 9b, presenting a general improvement of predictions for D eff .D data is the measured D eff and D pred is the corresponding predicted value from the updated Re TKE model.This newly fitted c provides an insight to correct the coherent structure contribution in the dense vegetation canopy (ah > 1.0) due to limited spaces for the structure to develop.However, as discussed in Tseng and Tinoco (2021), this empirical coefficient could vary under various flow conditions and environmental constraints.More experimental or numerical studies are required in the future to obtain a more comprehensive relation between c and ah.

Conclusions
We conducted a series of experiments using acrylic cylinders to construct different layouts of simulated vegetation canopies to study the effect of stem size and spatial heterogeneity on the interfacial gas transfer process across AWI and SWI in a meter-scale race-track, recirculating flume.The results were compared with data from previous experiments in flows with and without vegetation.Existing models for predicting interfacial gas transfer across AWI and SWI were updated by implementing a newly proposed method to quantify the randomness level of the vegetation array using the spectra of spatially distributed quantities.This study provides important insight to account for vegetation distribution patterns in laboratory-based empirical AWI and SWI gas transfer models to improve predictions for more universal applications to full-scale field conditions.Major contributions of this study are concluded as follows: • According to PIV data, mean flow and TKE production temporally and spatially averaged over a horizontal plane are representative for emergent vegetated flows, showing that the horizontal-shear stem-wake turbulence is dominant in the system.However, in submerged canopies, mean flow and TKE production spatially and temporally averaged over a vertical plane captures the representative peak value on top of the canopy, indicating that vertical-shear canopy-scale eddies are dominant for the transfer process.• We propose a newly defined randomness indicator, ξ, according to the corresponding FFT spectra of 〈u〉 2 and 〈u〉 3 .The approach decomposes the spanwise flow pattern into its wavenumber components, which provides the information about the periodicity characteristic and pattern regularity to quantify the randomness level of the canopy.• By using the horizontal-shear mean TKE production, P uv,m , as the predictor in the modified SR model to predict air-water surface gas transfer in emergent vegetated flows, stem size and canopy layout do not have significant influences on the model results.However, when using the maximum streamwise-averaged verticalshear TKE production, P uw,max , as the predictor for the submerged cases, we need to implement the newly proposed η correction parameter based on the randomness indicator, ξ(〈u〉 3 ) , to determine how representative the peak value is for the whole array.The significance of this correction arises from the method of measuring P uw,max based on the peak value occurring atop a specific row of the vegetation canopy, which staggered setup for laboratory simplifications could cause wrong estimates of gas transfer rates.• For the interfacial transfer process across SWI, the Re TKE -dependence model proposed by Tseng and Tinoco (2022) was modified to more accurately estimate the total near-bed TKE contribution from the temporal varying coherent structures by the non-linear stem-flow-bed interactions in random canopies.Using the proposed η correction parameter based on the randomness indicator, ξ(〈u〉 2 ) , the empirical coefficient, c, for fully developed coherent eddies was corrected.Furthermore, in the dense vegetation canopy (ah = 1.0), the limited space between vegetation elements restricts the development of the coherent vortex structures, where the viscous layer around the surrounding stems' surface would reduce the vortex strength and break the vortices into small isotropic eddies.Hence, c needs to be further corrected by a smaller fitted value (c = 0.9), to estimate a contribution from the coherent structures in the dense case for a more precise prediction.
The proposed updated interfacial transfer models provide a predictive framework for water quality management and ecosystem restoration, under a more realistic natural environment with varying stem sizes and canopy layouts in vegetated flows.To apply the model in a real-world scenario, it is necessary to obtain cross-stream velocity measurements on the horizontal plane with the same level of resolution as those measured in the corresponding staggered configuration set up in the laboratory (as described in Section 3.3).Consequently, the field measurements (random) and the laboratory measurements (staggered) can be directly compared, ensuring the acquisition of a meaningful ξ value for the η correction in the updated model (Equation 29).

Water Resources Research
10.1029/2023WR035220 So far, the approach has been proven to be applicable to fields in which highly resolved cross-stream slices of mean velocity can be measured and compared to the laboratory results, indicating opportunities for future research in exploring potential optimizations and broader applications to more complex scenarios.For instance, it is worth investigating more detailed influences from different layouts of the canopy.More cases with varying randomness of the vegetation distribution would help test and improve the current updated interfacial transfer models for broader range canopy configurations, while potentially developing an empirical way to parameterize ξ with relevant array characteristics and hydraulic parameters.Furthermore, a more robust scaling method can be developed to estimate the minimal space required for coherent structures to fully evolve, which could be related to the stem size, d, inter-stem spacing, bottom roughness, k s , and stem-scale Reynolds number, Re d .
Simplifications are applied in the experiments to gain a clearer understanding of the underlying principles and mechanisms, so as to update the existing gas transfer models and include the array randomness effects.Such simplifications indicate some unresolved problems, which provide further opportunities for future studies.For instance, the first-order reaction model for bio-chemical consumption (Equation 6) is a simplified assumption as the current study focused on only the long-term mean gas transfer rate.However, such relations between flow strength, bio-film growth, and bio-chemical consumption can be more complex and sensitive to instant gas transfer locally, which shows another future research direction.Also, experiments in the current study were set within a specific range of hydraulic and sediment parameters to avoid significant bedform changes that might affect the SWI gas transfer and complicate the system.Hence, future work on bedform dynamics and the corresponding effects of the near-bed flow turbulence and the SWI gas transfer under different vegetation layouts could be a great extension.CT acknowledges funding support from Taiwan-UIUC Fellowship during his Ph.D. study at UIUC.This study was supported by NSF through CAREER EAR 1753200.Any Opinions, findings and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect those of the National Science Foundation.The authors are grateful to the associate editor and two anonymous reviewers for their thoughtful and constructive comments on improving and clarifying the article.

Figure 1 .
Figure 1.(a) Schematic of the conceptual model of gas transfer across AWI and SWI in vegetated rivers.U is the mean streamwise flow velocity, δ c is the diffusive layer thickness (not to scale), C DO is the spatially averaged DO concentration away from the bed, C DO,b is the near-bed spatially averaged DO concentration, and C w is the local instantaneous DO concentration.The figure is adapted from Figure 1 in Tseng and Tinoco (2022).(b) Schematic of the plan view velocity transect (measured along the vertical dashed line) within a staggered (above the horizontal dashed line) and a random (below the horizontal dashed line) vegetation array.

Figure 2 .
Figure 2. (a) A 3D overview sketch of the race-track flume.(b) Side-view sketch of the straight test section (not to scale), where the water depth, H, ranges from 10 to 20 cm, the thickness of the sediment bed, H s = 10 cm, and the stem diameter varies from, d = 0.635 cm to d = 2.54 cm.The red dashed area indicates the observation window for PIV measurement and the blue dashed area shows the region used for calculating the x-direction-averaged values.Two DO sensors were put near the surface and near the bed at the upstream side of the vegetation array as shown by the red arrows.The inset shows the image of the near-bed sensor.(c) Partial region of the two vegetation array layouts (random and staggered).Two configurations were designed under the same stem density.The green dashed line shows the location of the vertical PIV laser slice, focusing on the center row of the array.(d) The race-track flume picture (right) and top-view experiment images of the thick (top-left: d = 2.54 cm) and thin (bottom-left: d = 0.635 cm) vegetation patches.

Figure 3 .
Figure 3. (a) The re-aeration curves obtained from the DO measurements near the bed (red dashed line) and near the surface (blue dashed line) under staggered vegetation array with d = 0.635 cm, h/H = 1.0,U m = 11.1 cm/s (refer to case No. 12 in Table 1).The bulk, near-bed DO concentration difference, ΔC DO , is defined as the difference between the spatial-averaged DO concentration away from the bed, C DO , and the spatial-averaged near-bed DO concentration, C DO,b , which is shown by the yellow shaded area.(b) The ΔC DO curve during the re-aeration process corresponding to (a).

Figure 4 .
Figure 4.The normalized time-averaged flow field, mean flow velocity profile, streamwise-averaged TKE profile, and streamwise-averaged P uw profile at the vertical plane under random and staggered array layouts for stem diameters d = 0.635 and 2.54 cm in submerged vegetation canopy (h/H = 0.5).The corresponding U and mean flow Reynolds number, Re H , for cases (a) to (d) are: U = {6.31;5.90; 4.81; 4.06} cm/s, and Re H = {8,303; 7,764; 6,630; 5,349}.The black dashed line represents the height of the vegetation canopy.The sketch shows the side-view and plan-view of the vertical PIV setup, where red dashed area shows the observation region for PIV measurement and the blue dashed area shows the region used for calculating the streamwise-averaged values.

Figure 5 .
Figure 5.The normalized time-averaged flow field and mean flow velocity profile at the vertical plane under random and staggered array layouts, with stem diameters d = 0.635 and 2.54 cm in emergent vegetation canopy (h/H = 1.0).The corresponding U and mean flow Reynolds number, Re H , for cases (a) to (d) are: U = {5.16;8.02; 6.27; 3.86} cm/s, and Re H = {4,273; 4,320; 4,689; 6,791}.The sketch shows the side-view and plan-view of the vertical PIV setup, where the red dashed area shows the observation region for PIV measurement and the blue dashed area shows the region used for calculating the streamwise-averaged values.

Figure 6 .
Figure 6.The normalized time-averaged flow field, mean flow velocity profile, TKE profile, and streamwise-averaged P uv profile at the half-stem-height horizontal plane (z = h/2) under random and staggered array layouts with stem diameters d = 0.635 and 2.54 cm in emergent vegetation canopy (h/H = 1.0).The corresponding horizontal-mean flow velocity, U m , and the stem-scale Reynolds number, Re d , for cases (a) to (d) are: U m = {5.38;5.44; 5.91; 8.55} cm/s, and Re d = {342; 346; 1,500; 2,173}.The sketch shows the side-view and plan-view of the vertical PIV setup, where the red dashed area shows the observation region for PIV measurement and the blue dashed area shows the region used for calculating the streamwise-averaged values.The above results only showed the center spanwise region (≈±20 cm from the centerline, y = 0), as the y-axis was cropped from y/d = 30 to 30 when d = 0.635 cm and y/d = 8 to 8 when d = 2.54 cm.

Figure 8 .
Figure 8.The experimental data from the current tests and from previous studies (Tseng & Tinoco, 2020, 2022) and the prediction of air-water surface gas transfer rate, k L , by the newly modified SR model in (a) emergent canopy and (b) submerged canopy.In the emergent case (a), the critical stem-scale Reynolds number, Re dc = 300, and the best-fitted experimental coefficient, α = 0.12, when Re d is above Re dc .The inset shows the zoomed-in region in the red box.In the submerged case (b), the critical upper-canopy Reynolds number, Re upc = 25, 000, and the best-fitted experimental coefficients, α = 0.065 and α = 0.71, when Re up is below and above Re upc .Blue markers represent data obtained under random distributed arrays, R. Red markers represent data obtained under staggered distributed arrays, S.
Figure9ashows the modified Re TKE -dependence model prediction byTseng and Tinoco (2022) for the SWI gas transfer diffusion coefficient, D eff , against previous observations and the current experimental data.With the corrected k tcs term based on Equation 37, data from cases with d = 2.54 cm (ah = 0.25) (black crosses) generally follows the previous experiments and the model prediction line.However, data from cases with d = 0.635 cm (ah = 1.0) (red circles) shows a slight deviation from the predicted line and other data points under sparse canopies (black crosses: d = 2.54 cm (ah = 0.25), blue circles: Tseng and Tinoco (2022) (ah = 0.5 and 0.1)), showing an overestimation of Re TKE when the canopy is dense.This is because, in the dense canopy (ah = 1.0), the limited space between vegetation elements restricts the development of such temporal-varying vortex structures, reducing the strength of these coherent systems closely surrounded by the viscous layer around the stems' surface.Under this condition, the originally downstream-stretching vortices would be broken down into smaller eddies(Pope, 2000) with a weaker component from the k tcs term, suggesting that the coefficient c would need to be updated with a smaller value.Replacing c with a new best-fitted value c = 0.9 for the dense canopy cases (d = 0.635 cm, ah = 1.0), all data from channels with vegetation under different stem sizes and array densities (color markers) now matches better with data from bare-bed channels (gray markers) and follows the updated model predictions as shown in Figure9b.The total normalized variation was calculated as:

Figure 9 .
Figure 9.The randomness-corrected Re TKE -dependence model prediction for the SWI gas transfer diffusion coefficient, D eff , (a) without and (b) with the spatial correction for the development of instantaneous coherent structures against previous observations and laboratory data from the current experiments.The three regimes: molecular, dispersion, and turbulence are separated by the dotted vertical lines.Data under vegetated flows are represented by the color markers as follows -red circles: current experimental data with d = 0.635 cm and ah = 1.0, black crosses: current experimental data with d = 2.54 cm and ah = 0.25, blue circles: Tseng and Tinoco (2022).Previous data under bare-bed open-channel flows are represented by the gray markers as follows -circles: Richardson and Parr (1988), stars: Marion et al. (2002), diamonds: Nagaoka and Ohgaki (1990), squares: Lai et al. (1994), crosses: Packman et al. (2004), triangles Chandler et al. (2016), and inverted triangles: Roy et al. (2004).The insets in (a, b) only show the experiment data under vegetated flows and the Re TKE model prediction line.V t is the value of total normalized variation calculated by Equation 38.

Table 1
Experimental Parameters: Array Configuration, Array; Stem Diameter, d; Solid Volume Fraction, ϕ; Roughness Density, ah; Submergence Ratio, h/H; Time-Averaged Mean Flow Velocity Along the Centerline, U; Time-Averaged Mean Flow Velocity Along the Centerline Within the Vegetation Canopy (Only Consider the Region From the Bed Up to the Vegetation Height), U veg ; Time-Averaged Horizontal-Mean Flow Velocity at z = 0.5 hr, U m ; Mean Flow Reynolds Number, Re H = (U m R H )/ν for Emergent Arrays and Re H = (UR H )/ν for Submerged Arrays, Where R H Is the Hydraulic Radius; The model distinctly shows three dominant regimes based on Re TKE :