Dune Geometry and the Associated Hydraulic Roughness in the Fluvial to Tidal Transition Zone of the Fraser River at Low River Flow

In deltas and estuaries throughout the world, a fluvial‐to‐tidal transition zone (FTTZ) exists where both the river discharge and the tidal motion drive the flow. It is unclear how dune characteristics are impacted by changes in tidal flow strength, and how this is reflected in the hydraulic roughness. To understand dune geometry and variability in the FTTZ and possible impacts on hydraulic roughness, we assess dune variability from multibeam bathymetric surveys, and we use a calibrated 2D hydrodynamic model (Delft3D‐FM) of a sand‐bedded lowland river (Fraser River, Canada). We focus on a period of low river discharge during which tidal impact is strong. We find that the fluvial‐tidal to tidal regime change is not directly reflected in dune height, but local patterns of increasing and decreasing dune height are present. The fluvial‐to‐tidal regime change is reflected in dune shape, where dunes have lower leeside angles and are more symmetrical in the tidal regime. The calibrated model allows to estimate local patterns of dune heights using tidally averaged values of bed shear stress. However, the spatially variable dune morphology hampers local dune height estimation. Changes in dune shape do not significantly impact the reach‐scale roughness, and estimated dune roughness using dune height and length is similar to the dune roughness inferred from model calibration. Hydraulic model performance with a calibrated, constant roughness is not improved by implementing dune‐derived bed roughness. Instead, the data analysis revealed that large‐scale river morphology may explain differences in model roughness and corresponding estimates from dune predictors.

• Hydraulic roughness in the fluvial-to-tidal transition zone estimated from dune geometry agrees with calibrated model roughness • Variation in dune asymmetry and leeside angle across a fluvial-to-tidal transition zone has little impact on reach-scale hydraulic roughness • Estimated spatial bedform patterns from modeled shear stress match measured bedform patterns, but absolute dune heights do not

Supporting Information:
Supporting Information may be found in the online version of this article. 10.1029/2023JF007340 2 of 25 how bedforms and their corresponding roughness respond to a change in hydraulic regime, even though dune geometry and roughness prediction are essential for river management (ASCE Task Force, 2002;Best, 2005;Warmink et al., 2013), interpreting sedimentary rock records (Das et al., 2022), and understanding sediment fluxes (Venditti & Bradley, 2022).
Dunes are rhythmic undulations of the river bed with heights of approximately 10 −1 -10 0 m and lengths around 10 1 -10 2 m, thereby being larger than ripples but smaller than bars.They adjust to changes in the hydraulic regime, but not in a consistent manner.Until recently, it was often assumed that any spatial variability in dune geometry was caused by dunes that are not in equilibrium (Bridge, 2003;Carling et al., 2000;Holmes & Garcia, 2008), and the primary geometry (dune height and length) of equilibrium dunes was assumed to scale with water depth (Yalin, 1964).However, recent research has shown significant local spatial variation in dune height (Bradley & Venditti, 2017;Murphy, 2023;Venditti & Bradley, 2022) in riverine systems, including the Fraser River (Murphy, 2023), independent of the water depth.In the FTTZ, this variability is expected to be even more pronounced, since tidally influenced currents impose spatially varying water level fluctuations (and therefore bed shear stress changes) on diurnal and semi-diurnal time scales (Hoitink et al., 2003;Sassi et al., 2011).The resulting spatial longitudinal variability of dune geometry in the FTTZ is understudied, but two key studies exist.Prokocki et al. (2022) studied dunes in the lower 90 km of the Lower Columbia River (USA), and recognized differences in shape and 3D planform of dune geometry across the study reach.They related the changing dune morphology to downstream variations in grain size and spatiotemporal changes in tidal and fluvial flow.In the thalweg, they observed small upstream-oriented dunes downstream, and larger downstream-oriented dunes upstream.Unfortunately, they did not report on flow conditions in those distinct regions, or on the resulting hydraulic roughness differences.Lefebvre et al. (2021) studied 4-year long bathymetric data of the downstream 160 km of the Weser Estuary in Germany.They did not observe a clear trend in dune geometry in the longitudinal direction, but found dunes that are generally smaller than predicted based on the water depth (Allen, 1978).They did not provide information on the flow conditions or resulting roughness.Beyond these recent studies, the response of dune geometry to shear-stress variation at the change from a fluvial to tidal regime is unknown, and it is uncertain if dune geometry predictors apply here.
To date, it remains unclear to what extent variability in dune geometry is relevant for the large-scale roughness needed to model hydrodynamics in the FTTZ.Bedforms, especially dunes, are known to be a major source of roughness in lowland rivers (Gates & Al-Zahrani, 1996;Julien et al., 2002).Hydraulic roughness quantifies the resistance to flow by objects protruding into the water column (Chow, 1959).It consists of skin friction drag (grain resistance) and form drag, resulting from pressure variation (Einstein, 1950).Dune variability can impact roughness parameterizations via its impact on the form drag.When modeling the FTTZ hydraulically, a roughness value must be chosen, generally expressed as Manning's n, the Chézy coefficient, or the roughness height.We identify three options to adopt at a suitable roughness value.First, the apparent hydraulic roughness can directly be inferred from calibration.The calibrated roughness is often represented by a single constant coefficient (Paarlberg et al., 2010); limited examples are available where a variable roughness coefficient is used.De Brye et al. (2011) used a linearly decreasing roughness coefficient in a model domain from a delta apex to the coast, to represent a gradual transition from the riverine to the marine environment.However, proof of the validity of this approach is lacking.
Second, a roughness predictor, such as that of Van Rijn (1984), can be used to estimate hydraulic roughness, which have been developed mainly in alluvial conditions without tides.With an increasing amount of bathymetric data available from field surveys, measured dune geometry can be used as input for these predictors.Typically, even in this option for hydrodynamic modeling, a calibration constant is still implemented to adjust the estimated hydraulic roughness.Third, when no field data is available, such as for inaccessible regions or during peak flow conditions, roughness predictors can be employed using dune characteristics estimated from predictors such as Van Rijn (1984), Yalin (1964), Karim (1995), and Venditti and Bradley (2022).The latter approach needs iteration, since roughness impacts the water levels, which in turn impacts the dune geometry estimations.The reliability of these dune geometry predictors in the FTTZ is unknown, and needs testing.There is a need to improve hydraulic roughness parameterization and dune geometry estimations in the FTTZ, since the output of hydrodynamic models strongly depends on an accurate specification of roughness (Lesser et al., 2004;Morvan et al., 2008;Wright & Crosato, 2011).
In this research, we aim to increase understanding of bedform variability and related roughness that occurs at the transition from a shallow mixed tidal-fluvial regime to a tidal regime.To do so, we assess the bedform 10.1029/2023JF007340 3 of 25 characteristics and the resulting hydraulic roughness of the FTTZ of the Lower Fraser River.The Lower Fraser River is a sand-bedded lowland river in British Columbia, Canada, with a significant decrease in tidal energy (amplitude) 40 km landward of the river mouth (Wu et al., 2022).We aim to answer three research questions.

Field Site
The Fraser River (Figure 1) is located in British Columbia, Canada, and drains 228,000 km 2 of mountainous terrain.
The Fraser exits a series of bedrock canyons approximately 185 km upstream of the river mouth at river kilometer (RK) 0, where it enters the gently sloping Fraser Valley, past the towns of Hope (RK 165) and Mission (RK 85).
The Fraser River has an annual river discharge of 3,410 m 3 s −1 at Mission (Water Survey of Canada (WSC) Station 08MH024), but flow rates vary between a mean low discharge of 1,000 m 3 s −1 in winter time (November-April) and a mean high discharge of 9,790 m 3 s −1 during the snow melt-dominated freshet in May, June and early July (Attard et al., 2014;McLean et al., 1999).At New Westminster, 34 km upstream from the river mouth, the river bifurcates, forming the Fraser Delta where the Main Channel splits into two tributaries: the North Arm and the Main Channel.
The Main Channel carries 88% of the flow, until it bifurcates into Canoe Pass (RK 13), which conveys approximately 18% of the total flow (Dashtgard et al., 2012;NHC, 2008;WCHL, 1977).The fluvial-to-tidal transition zone of the river is influenced by a predominantly semi-diurnal tide (Wu et al., 2022), with a mean tidal range at the mouth of approximately 3 m (Kostaschuk & Atwood, 1990).The tidal motion influences water levels up to Mission during high flow, but can penetrate up to Chilliwack (RK 120) during base flow conditions in August till April, creating a strong backwater effect (Kostaschuk & Atwood, 1990) (i.e., an increase in subtidal water level as response to the tidal forcing).The Fraser has an undammed, unregulated flow, which reflects climatic conditions.Human-made influences include dikes (90% of the reach), pipelines and bridge constructions, and dredging of the Main Channel occurs.
The Port of Vancouver dredges from the river mouth (RK 0) to the Port Mann Pumping Station (RK 42), with the most significant dredging in the deltaic reach from RK 35 to the river mouth (Nelson, 2017) to maintain a constant fairway depth (McLean & Tassone, 1989;Stewart & Tassone, 1989).The depth is larger in the tidal region, and has been made deeper by dredging.This results in a large decrease in momentum flux (Wu et al., 2022), defined as water mass times velocity in a unit area, at the upstream limit of the dredging works.Additionally, Wu et al. (2022) related this decrease in momentum flux to the influence of the Pitt River.They identified the junction of the Pitt River as the transition from a tidally dominated to a river-dominated regime, and they noted the importance of this system for tidal attenuation.Therefore, two different regimes can be identified in the study area.The first regime, hereafter called the tidally dominated regime, occurs seaward of RK 40, and is characterized by a strong influence of tides and a large tidally averaged water depth compared to the more landward-located region.The second regime is the mixed tidal-fluvial regime, in which tides are less strong and the water depth is shallower, which occurs landward of RK 40.This roughly coincides with the upstream end of the modern Fraser Delta at RK 35 (Venditti & Church, 2014;Venditti et al., 2015), where the Fraser River bifurcates into the North Arm and the Main Channel.
The difference in tidal strength in the two regions does not impact grain sizes of bed sediments in the thalweg.
The main transition of grain size characteristics occurs around RK 100.Upstream of RK 100, the bed of the Fraser River consists of gravel, and downstream of a gravel-sand-transition near Mission, the main bed material is sand (Venditti & Church, 2014).The median grain size (D 50 ) is 351 μm and the mean is 415 μm (Figure 2).There is a minor trend of downstream fining in the thalweg of the lower 50 km of the river (1.14 μm per kilometer, Figure 5c), resulting on average in a decrease in D 50 of approximately 100 μm over this reach, although there is considerable scatter due to gravel and mud deposits along the banks.The data underlying this figure is a compilation of multiple sources.The samples up until RK 48.5 were collected by McLaren and Ren (1995), who sampled bed material in the Main Channel and delta front at 0.5 km increments with a Shipek sampler.Although this grain size data is decades old, broad patterns are likely to be consistent with present conditions (Venditti & Church, 2014), and grain size shows little seasonal or year-to-year variation (Kostaschuk et al., 1989;McLean et al., 1999;Pretious, 1956).Venditti and Church (2014) measured 33 samples of RK 48.5-80, with a dredge sample in 2009, and Murphy (2023) collected 115 additional samples in this same reach using a pipe dredge.
They did not perform analysis on the fraction smaller than 63 μm.The Pitt system does not impact the sediment composition of the Fraser, since the net bedload transport is directed upstream toward Pitt lake (Ashley, 1980).In the delta, the river deposits its sand load in the channel and its banks, and its silt load on the distal margins and tidal flats (Venditti & Church, 2014) (Figures 2a-2c, and 2d).Seaward of the river mouth, the grain size decreases dramatically to a D 50 of 22 μm.Locally, the river interacts with its bank and bed substrate.Gravel and clay patches are present at the outer banks on the north and south sides of the river.These patches are either modern deposits, such as gravel bars, or older Pleistocene glacial deposits seen as fine outwash deposits and coarse glacial till (Nelson (2017), Figure S1 in Supporting Information S1), constraining the river's course.
This study focuses on the Main Channel of the Fraser River, from the confined part of delta mouth at Steveston Harbor at RK 10, to Mission at RK 85 (Figure 1).The area is located in the FTTZ, and low-angled dunes (Bradley et al., 2013;Kostaschuk & Best, 2005;Kostaschuk & Villard, 1996), with no or intermittent flow separation, cover the river bed.When assessing local scale processes, we focus on three zones in the FTTZ (Figure 2; Figure S7 in Supporting Information S1).The zones are located at RK 21.5-23 (zone 1; tidal regime), 50-52.5 (zone 2; fluvial-tidal regime) and 57-59.5 (zone 3; fluvial-tidal regime).The selection of zones is based on three criteria.First, a decreasing amount of tidal energy (tidal amplitude) from zone 1 to 3. Second, a simple geometry, without bifurcations or confluences, to limit the influence of complex currents on dune geometry.Third, a limited degree of human influence on the river bed.Zone 1 is 1 km shorter than the other zones due to dredging along the downstream side, and an engineering structure on the upstream side.Markers differentiate between samples taken in the thalweg (solid marker) or outside along the river banks (indicated with "x").The data is from a data compilation by Venditti and Church (2014), which includes reanalyzed observations from McLaren and Ren (1995), and recent observations by Murphy (2023).

Field Data and Preprocessing
Raw multibeam echosounder (MBES) riverbed data were provided by the Public Works and Government Services, Canada.The measured bathymetry comprises data of the Main Channel between river kilometer −1 to 85 and covers the navigation channel, but does not provide full bank-to-bank spatial coverage.Data were collected during base flow conditions in January, February and March of 2021.This period is characterized by relatively constant discharge and little change to the bed surface (Bradley & Venditti, 2021).During the survey period, the measured discharge (at an hourly interval) was relatively constant, with monthly mean values of 1,416 m 3 s −1 (standard deviation (SD) 184 m 3 s −1 ), 1,051 m 3 s −1 (SD 140 m 3 s −1 ) and 1,074 m 3 s −1 (SD 35 m 3 s −1 ) at Hope (RK 165) for the 3 months, respectively (Water Survey of Canada, Station 08MF005).
The MBES data is gridded onto a 1 × 1 m 2 grid, and the resulting MBES data sets contain x, y, and z coordinates.
Next, all bed level data is converted from Cartesian (x, y) coordinates to curvilinear coordinates (s, n) with the same spatial resolution (Vermeulen, Sassi, & Hoitink, 2014).Herein, s is the streamwise direction, and corresponds with river kilometers (RK, i.e. chainage) and n is the stream normal direction, wherein n = 0 m is defined as the central river axis, which roughly coincides with the thalweg.

Hydraulic Model Setup
A 2DH (two dimensional horizontal) hydraulic model was set up in the Delft3D Flexible Mesh (FM) model suite (Kernkamp et al., 2011).The model simulates depth-averaged flow quantities based on the two-dimensional shallow water equations.The numerical domain covers the Fraser River from river kilometer 85, to the offshore region of the Strait of Georgia reaches where depth exceeds 200 m.The bathymetry of the model of Wu et al. (2022) was taken as a basis, and the higher resolution MBES data described above were used for the bathymetry of the channels in the estuary.The MBES data was smoothed with a LOESS filter (Cleveland, 1979;Cleveland & Devlin, 1988;Ruijsscher et al., 2018), and interpolated on an unstructured curvilinear grid.The median cell size was 30 m, and varies between 5 m in the river and 1,000 m offshore.Bars that do not get submerged during an average yearly freshet (flood) were not well-represented in the bathymetric data, so bar elevation in the model was manually increased to 10 m above mean sea level, to prevent flooding.
The upstream boundary condition is the discharge at Mission (RK 85) for the time period of November 2017 to October 2018.The discharge at Mission is estimated using a rating curve if the discharge exceeds 5,000 m 3 s −1 .At lower discharge conditions, tidal influences make the rating curve at Mission inaccurate, and therefore the discharge was calculated using the discharge at Hope (RK 165) and two smaller tributaries (Chilliwack River and Harrison River).Using this calculation method, discharges measured at Mission (larger than 5,000 m 3 s −1 ), were on average underestimated by 3%, and no significant temporal delay was observed.Additionally, at two downstream confluences, a constant discharge of 315 m 3 s −1 at Stave River (RK 74) and 130 m 3 s −1 from the Pitt River (RK 45) were added to the Fraser flow.Stave River is dammed at 3 km upstream, therefore having a controlled flow.The Pitt River drains a lake and has therefore a nearly constant outflow.At the downstream boundary (see Figure 1b), water levels influenced by tides are imposed.Eight primary tidal constituents, the most important overtide (M4) and compound tides are determined via the Delft Dashboard toolbox (Van Ormondt et al., 2020), using the TPXO8.0database (Egbert & Erofeeva, 2002).
The tidal amplitudes and phases at the downstream boundary were corrected to minimize the error in the modeldata comparison at the Point Atkinson tidal gauging station.This correction was on average 0.8% of the tidal amplitude and 20° for the tidal phases, for the 13 tidal components.The model was calibrated for low discharge (<1,600 m 3 s −1 ; Figure 3b), to simulate flow conditions that correspond to the low-discharge bathymetry.The calibration was performed by varying the Manning's roughness coefficients and evaluating the resulting water levels and tidal amplitudes of the three most important tidal constituents at 7 gauging stations (RK 0,9,18.5,35,42,70,85)(Figure 2).The principal tidal constituent M 2 is used for calibration, together with M 4 and K 1 .Relative phase differences between M 2 and M 4 (the first overtide of M 2 ) influence tidal duration asymmetry, the main mechanism for driving residual bed-load transport in estuaries (Van De Kreeke & Robaczewska, 1993).The diurnal tide K 1 is relatively large at the west coast of North America, and interaction between diurnal and semi-diurnal frequencies can produce asymmetric tides as well (Hoitink et al., 2003).
The tidal amplitudes were derived from harmonic analysis using t_tide (Pawlowicz et al., 2002).The best performing model had a uniform Manning's roughness coefficient of 0.026 s m −1/3 (Figure 3).Disconnecting the hydraulic roughness at the regime transition at RK 40, thereby allowing for two different roughness values, did not improve the calibration (Figure S3a in Supporting Information S1).The uniform Manning's coefficient (roughness) of 0.026 s m −1/3 is slightly higher than the calibrated Manning's coefficient of Wu et al. (2022), who used a uniform roughness of 0.015 s m −1/3 .The difference in roughness value is attributed to the difference in grid resolution.Our model grid in the river domain is overall coarser than the model of Wu et al. (2022) who used a 10 m resolution, so that the schematization of the bathymetric data on our grid results into slightly wider channels.Our value for roughness is considered to be more appropriate for natural sand-bedded rivers (Chow, 1959).

Analysis of River Bathymetry: Dune Detection
Bathymetry was analyzed to derive dune characteristics.Three longitudinal profiles were taken, along the centerline and at approximately 80 m from the north and south bank, and their longitudinal location is expressed as chainage (km) that is, along-river kilometer.In three focus areas (Figure 1), a longitudinal profile was taken every 10 m, resulting in 27, 41, and 23 profiles for zones 1, 2, and 3 respectively, depending on the river width.
To ensure that bedforms in all profiles were primarily caused by natural mobile bed conditions, we excluded bathymetry that showed extensive scour, human-made structures and dredging marks (Figure S2 in Supporting Information S1).
From the filtered profiles, bedform characteristics were determined by using a standard bedform tracking tool (van der Mark et al., 2008).In the tool, the default filter span (c = 1/6) was suitable to filter out small features such as measurement errors or outliers.Three span values (P0), corresponding with bedforms with a length of 20 m ± 10, 50 m ± 20 and 100 m ± 30, were used as input to detrend and smooth the profile.The span values in the tool are based on a spectral analysis yielding the dominant bedform wave lengths in each section.
Based on the detrended and smoothed profile, a zero-crossings profile was defined, based on which individual dunes were identified, and dune characteristics were calculated (see van der Mark et al. ( 2008) for more 8 of 25 information).Dune characteristics include dune height Δ (m), the vertical distance between the crest and downstream trough, dune length λ (m), the horizontal distance between two subsequent crests, the seaward-facing leeside angle LSA (°), the slope from a linear fit of the dune's leeside, excluding the upper and lower 1/6 of the dune height, and the landward-facing stoss side angle SSA (°) calculated in the same manner as the leeside angle.
The bedform slip face angle SFA (°), the steepest part of the leeside, and is defined as the 95-percentile of the leeside angle.Finally, bedform asymmetry is calculated as the ratio between the length of the (seaward) leeside and the total bedform length (Lefebvre et al., 2021).
Bedforms with heights smaller than 0.1 m are not distinguishable from the error of the survey, and are excluded from the analysis.Bedform lengths smaller than 3 m are excluded as well, since the resolution of the bathymetric data (1 m) is too small to detect small bedform features.Features with a height greater than 2.5 m (2% of all detected bedforms) or a length greater than 200 m (0.08% of all detected bedforms) are classified as bed features unrelated to dunes such as scour holes or wake deposits downstream of human-made structures.These had a different geometry than mobile dunes, which was confirmed by visual inspection of the bathymetric data.Large dunes (>2.5 m) as found in previous studies (Kostaschuk & Luternauer, 1989;Pretious & Blench, 1951;Venditti et al., 2019) were not observed in the low-discharge bathymetry used in this study.

Analysis of River Geometry
To assess if and how river geometry impacts dune height and length, river geometry was parameterized by river width, curvature, transverse bed slope and excess depth.River width W (m) was determined from a polygon following the longitudinal low water line, which is considered to be the discharge carrying section of the river under low discharge conditions.Cross-sectional area A (m 2 ) was subsequently approximated from the tidally averaged water depth and river width, assuming a trapezoidal shape of the cross-sectional area, where the river bank has a 60° slope.Curvature r (km −1 ) was defined as the inverse of bend radius following the approach of de Ruijsscher et al. (2020).Local transverse bed slope ξ (−) was defined as the slope between the two sides of the main river channel, longitudinally discretized at intervals of 100 m.Finally, an excess depth parameter D e (−) was used as a measure to identify pools and scour holes (Vermeulen, Hoitink, et al., 2014), defined as: where D r is the regional mean depth of a discretized section of 500 m long, and D max the local maximum depth in this section.Sign indicates the signum function, which returns the sign of a real number.

Analysis of River Hydrodynamics
To assess local flow conditions and tidal attenuation, the hydrodynamic model was evaluated during low flow conditions in March 2018 (Figure 3b).The flow magnitude and direction, water depth and bed shear stress magnitude and direction per grid cell were saved every 10 minutes in the simulation.All output data were tidally averaged using a Godin filter (Godin, 1972).The Godin filter removes the tidal and higher frequency variance to obtain a low-passed signal primarily caused by the river flow.
Besides transforming the data into along and across-channel coordinates (s,n) according to Vermeulen, Sassi, and Hoitink (2014), the flow and shear stress vectors were rotated, to transform their orientation to along-channel direction (s-direction).This allowed differentiation between the in-and outgoing currents, which are in s-and s-direction, respectively.The percentage of time that the flow reverses and flows upstream (reversal time t rev (%) was then calculated.

Dune Geometry Estimation
Flow data from the model was used to estimate dune height Δ (m), using several dune height predictors that include some measure of flow strength (Karim, 1995;Van Rijn, 1984;Venditti & Bradley, 2022;Yalin, 1964).Multiple predictors were used to assess their applicability in the FTTZ, since none of the predictors was explicitly developed for the FTTZ.Simultaneously, this allows assessment of potential differences of dunes in the FTTZ and dunes in the river.First, dune height was estimated using the widely applied dune geometry predictor of Van Rijn (1984).This predictor is based on 84 laboratory experiments, with grain sizes ranging from 190 to 2300 μm, and 22 field data sets (490-3,600 μm) of relatively wide rivers (width/depth >0.3) with unidirectional flow.This corresponds well with conditions found in the Fraser River, except that the Fraser experiences bidirectional currents.To account for this, values of water height and shear stress are tidally averaged, since bed-material sediment transport in the Lower Fraser River generally follows the pattern of mean velocity over the tidal cycle (Kostaschuk & Best, 2005).The tidal averaging is described in Section 3.3.3.Dune height is thus: in which D 50 is median grain size (m), h is the water depth (m) and transport stage T is: where u* is the shear velocity (m s −1 ), and   *  is the critical shear velocity (m s −1 ) for initiation of sediment movement.Shear stress (τ, N m −2 ) relates to shear velocity and is here expressed non-dimensionally as the Shields number (θ) as in: (5) in which g is the gravitational acceleration (9.81 m s −2 ), ρ s is the sediment density (2,650 kg m −3 for quartz) and ρ w is the water density (1,000 kg m −3 for fresh water).Therefore, Equation 3 can be rewritten as: When 50 μm < D 50 < 5,000 μm, the critical shields number θ c (−) can be approximated with a simple expression (Zanke, 2003).The resulting values of θ c are approximately 0.03 (medium sand): in which the particle Reynolds number Re p is: where the relative submerged density R = (ρ s − ρ w )/ρ w (−) and ν is the kinematic viscosity (m 2 s −1 ), which is slightly dependent on water temperature as ν = 4 * 10 −5 /(20 + t) in which t = temperature (°C).Here, ν = 1.3 * 10 −6 is used for 10°C.
We also estimate dune height using the predictor of Yalin (1964): The predictor of Karim (1995) builds on that of Van Rijn (1984) and Allen (1978).The predictor of Allen (1978) is not included in this research, since it is mostly based on laboratory experiments.The predictor of Karim (1995) is based on the suspension criterion which utilizes the shear velocity and the particle fall velocity (w s ): where w s follows from Ferguson and Church (2004) as: 11) in which C 1 and C 2 are constants with values of 18 and 1.0, respectively, for slightly irregular particles.

Hydraulic Roughness Estimation
To estimate the impact of dunes on the water levels in the study reach, the hydraulic roughness was determined.
The hydraulic roughness estimates were then compared to the calibrated model roughness, and a simplified version was implemented in the model.The total hydraulic roughness, expressed as a friction factor  f , results from form friction and grain friction (Einstein, 1950).Assuming dunes are the dominant structures causing form resistance, the total hydraulic roughness was estimated as in Van Rijn (1984): Herein, k s consists of form roughness height k sf and grain roughness height k sg : where D 90 is the 90th percentile of the grain size distribution, and where the calibration constant γ d is taken as 0.7 in field conditions (Van Rijn, 1984).Δ and λ are derived from field data.
In the modeling suite of Delft3D, roughness values of Manning's n, n man (s m −1/3 ), are converted to a Chézy coefficient C (m 1/2 s −1 ) via (Manning, 1891): in which R h is the hydraulic radius, which can be simplified to the water depth h (m) for rivers that satisfy W ≫ h.
The Chézy coefficient is converted to the dimensionless Darcy-Weisbach friction factor f m according to Silberman et al. (1963):

Hydraulics and Morphology of the Fluvial-To-Tidal Transition Zone
The tidally averaged water depth in the study area fluctuates between 3 and 18 m (Figure 4a).In the mixed-fluvial tidal regime of the river (RK > 40), it increases gradually in the seaward direction, and in the tidal regime (RK < 40) it remains constant.The local increase in water depth is reflected in the tidally averaged and instantaneous shear stress profiles (Figure 4b).The downstream-directed maximum shear stress increases from 0.4 N m −2 in the upstream area to 10 N m −2 at the river mouth.Similarly, the upstream-directed maximum shear stress in relation to flow reversal (indicated by a minus sign in Figure 4b) increases to 6 N m −2 .In contrast, the tidally averaged shear stresses remain relatively constant over distance (fitting a linear model gives a slope of 10 −5 N m −2 km −1 ).The tidally averaged shear stress is on average 0.64 N m −2 and fluctuates between −1.0 N m −2 (indicating an upstream directed shear stress at the most downstream area, RK 0) and 2.2 N m −2 (at RK 67).
The tidal effect on the water levels and flow direction weakens in the upstream direction, and the amplitude of the tidal constituents M 2 and K 1 decreases as the tides attenuate (Figure 4d).The M 2 tidal constituent shows a particularly strong decrease from RK 10 to 40, while landwards the tidal attenuation is minimal.In the most upstream reach, bidirectional currents can still be observed (Figure 4c).During low flow conditions, upstream (flood) flow occurs for 45% of the time at the river mouth (around RK 10), and decreases to 25% at the most upstream location of the study reach.
The morphology of the Fraser River does not show consistent trends in the streamwise direction.The river width fluctuates between 500 and 1,100 m (Figure 5a).The cross-sectional area of the river remains relatively constant in the more upstream part of the river (RK > 40), since river depth varies inversely with river width.The more downstream part features larger fluctuations in cross-sectional area, since water depth remains relatively constant (Figure 5a).The bed level (Figure 5b) shows large variations, but remains relatively constant in the downstream part.River curvature, transverse bed slope and depth excess are strongly related (R 2 = 0.15-0.61,p < 0.005, Figure 5c), which reflects the low-sinuosity meandering morphology.

Response of Dune Geometry to Tidal Hydraulics
Dune geometry in the study reach varies considerably (Figure 6), with dune heights up to 2.4 m (mean: 0.46 m, median: 0.39 m, SD: 0.28 m) and dune lengths up to 194 m (mean: 24 m, median: 16, SD: 22 m).Multiple scales of superimposed bedforms co-exist, although most of the bed is covered only by primary dunes.Patterns in dune geometry are apparent, with some areas of relatively low and short dunes, and others with increasing or decreasing dune heights and lengths.For example, the thalweg has large dunes around RK 68, 77 and 85, with increasing and decreasing dune heights upstream and downstream of those local maxima.Such patterns are not consistent over the whole river width however, and where relatively large dunes prevail on one part of the river (e.g., north side), dunes may be small on the opposite side (see for example around RK 68).This variation in dune height and length, along the cross-section and longitudinally, is expressed as the standard deviation of all dunes present in one unit of channel width.This allows for comparison between longitudinal and cross-sectional variability.The mean standard deviation in dune height and dune length in the cross-sectional direction (0.20 and 13.0 m, respectively) is twice as high as the variation along the longitudinal direction (on average 0.11 and 6.8 m, respectively) (Figure 6).
Local patterns in dune height and length are difficult to explain, and do not reflect the regime change around RK 40, trends in river geometry, or grain size in the thalweg.However, visual inspection reveals dune occurrence is primarily determined by grain size along the outer banks-when the grain size is too large (gravel) or too small (clay), dunes will be absent (Figure S2e and S2f in Supporting Information S1).When the river cuts into a clay or gravel layer on the north or south sides of the channel, abrupt changes in dune geometry can result.
The gradual increase in strength and duration of tidal currents in the seaward direction influences dune shape.First, the dune crests become sharper (Figures 1i-1k).Second, the leeside and slip face angle of the dunes decrease in downstream direction (Figures 7a and 7b).In particular the slip face angle decreases faster in the tidal regime than in the fluvial-tidal regime.Since the stoss side angle remains relatively constant (with a slight increase in the tidal regime), dunes become more symmetrical in the seaward direction (Figure 7c).When the asymmetry is equal to 0.5, dunes are perfectly symmetric.Dune symmetry at nearly every location up to a distance of 75 km from the river mouth, but becomes consistent at around 40 km from the river mouth, indicating the impact of the regime change.Results from a two-paired student t-test shows that the leeside angle, slip face angle and asymmetry is significantly different (at a 95% confidence level) in the tidal and the fluvial-tidal regime, while stoss side angle is not.The leeside angle directly correlates with flow-reversal time (Figure 7d) and maximum shear stress (Figure 7e), showing lower leeside angles and more symmetric dunes in seaward direction, although large variation is observed.

Dune Geometry Estimation From Model Output
Dune height predictors were applied to the FTTZ of the Fraser River at both reach scale and zone scale.The predictors were not specifically developed for tidal rivers with bidirectional currents, so input values were tidally averaged.
The predictor of Van Rijn (1984) works well when all data is reach-averaged (estimated Δ vRijn = 0.52 m, compared to measured Δ = 0.50 m; Figure 8a).However, it underestimates dune height in the mixed fluvial-tidal regime (by about 20 cm at RK 80), and overestimates it in the tidally dominated regime (by about 24 cm at RK 10).All other predictors are inaccurate and overestimate the dune height significantly, with an increasing error in the downstream direction (Figures 8b-8d).For example, the reach-averaged estimated dune heights are 0.87 m, 1.27 and 1.83 m for the predictors of Yalin (1964), Karim (1995), and Venditti and Bradley (2022), respectively.
Local variability in dune height in the study area is not captured in dune geometry predictors because of the considerable spatial variability in the measured dune geometry.To establish the degree to which local variability in dune properties relates to flow properties obtained with the 2DH hydraulic model assum- ing a constant roughness, we focus on three zones in the FTTZ (Figure 2).In those zones, flow characteristics are modeled (see Figures S7-S11 in Supporting Information S1) and the dune height predictor of Van Rijn (1984) (Equation 2) is applied to each zone using model output per grid cell.The dune height predictor of Van Rijn (1984) performs reasonably well in estimating the local spatial pattern of dune height in the three zones (Figures 9a-9c), but the mean dune height is overestimated for zone 1 and 3, and underestimated for zone 2. To assess the performance of the van Rijn model in predicting dune patterns, a bias correction is performed.Numerical values were added to or subtracted from the predicted dune height in order to minimize the RMSE, and assess the overall patterns in the dune field rather than the actual value (Figure S5 in Supporting Information S1).The bias-corrected RMSE values of dune height average 0.13 m, which indicates that the spatial pattern of dune heights is relatively well captured by the predictor.The dune height predictor of Van Rijn (1984) captures the main processes that determine dune height in tidal environments, but does not reliably predict absolute values.The dune height predictors of Yalin (1964), Karim (1995) and Venditti and Bradley (2022) are even worse on the local scale (Figures 9d-9l).Notably, the bias correction improves their performance (Figure S6 in Supporting Information S1), but they capture the pattern less than well than the predictor of Van Rijn (1984).

Comparison of Observed Dune Roughness and Model Roughness
Hydraulic roughness generated by dunes was calculated using Equation 13, which includes dune height and length, but does not include the leeside or stoss side angles.The estimated roughness decreases in the downstream direction (Figure 10), which is mainly caused by an increase in water depth.The main variability in roughness is due to variability in water depth, which is most pronounced in the upstream part (RK > 40) of the river (Figure 4).Additionally, local fluctuations in roughness correspond to the local patches of higher dunes, for example, at RK ∼ 54, 63, and 68.The decrease in grain roughness due to a subtle degree of downstream fining has a small impact on overall roughness, because grain roughness values are only around 3% of typical form roughness values.In the downstream reach (RK < 40), smoothed roughness shows a persistent out-of-phase relation with the gradient in smoothed bed elevation (moving average filter of 8 km) (Figure 10a) that is absent in the upstream part of the river.
The calculated dune roughness differs slightly from the uniform roughness used in the model (Manning's roughness of n man = 0.026 s m −1/3 ).The derived friction coefficient f m from the model's roughness (Equations 17 and 18) displays the expected decrease in seaward direction, reflecting the increase in water depth.Dune roughness agrees reasonably well to the uniform model roughness (Figure 10b), but local fluctuations are not well-represented.In  1964), (c) Karim (1995), and (d) Venditti and Bradley (2022).The measured data is based on the average of three longitudinal transects, and includes the minimum and maximum values in a gray shading.The modeled data is based on the model with a constant roughness of n man = 0.026 s m −1/3 .In subfigure (a) estimated dune height with dune-adjusted roughness (varying between 0.024 and 0.028 s m −1/3 ) is displayed in green.
the upstream region (RK > 40), the model roughness is slightly lower than the dune roughness, while they are similar in the downstream region.
To represent the effect of dune height variation on roughness in the hydrodynamic model, and to investigate if this can improve the calibration, the dune roughness as calculated by Equation 13is divided into three linear components: a uniform roughness of n man = 0.024 s m −1/3 from RK 0-10, a linearly increasing roughness of n man = 0.026-0.027s m −1/3 in the tidally dominated regime (RK 10-40), and a uniform roughness of n man = 0.028 s m −1/3 in the mixed fluvial-tidal regime (RK > 40).A small transition area between the breaks is implemented to ensure a smooth transition to the new roughness regime.These roughness transitions correspond with the transition from the fluvial regime to the deltaic regime around RK 40, and the downstream change from a confined to a less confined channel at around RK 10 (Figure 1).
The dune-derived roughness has little impact on the calibration of water levels and tidal amplitude of the M 2 , K 1 , and M 4 tidal components (Figure S3b in Supporting Information S1).On average, the RMSE value of the modeled water height decreases from 0.36 to 0.35 m, and the difference between the modeled and observed M 2 amplitude increases from 3% to 4% and K 1 from 4% to 6%.Additionally, using the dune-derived roughness for dune height estimation with the Van Rijn predictor only slightly improves the estimated values (17 cm at RK 80 and RK 10) (Figure 8).

How Are Bedform Characteristics Impacted by the Sudden Change in Tidal Flow Strength During Periods of Low River Discharge?
During the low river discharge conditions like studied in this research, the increase in water depth around RK 40 results into two different hydrodynamic regimes (Figure 11).The tidally dominated regime is characterized by large maximum absolute shear stresses, a large tidally averaged water depth, relatively symmetrical dunes, low leeside and slip face angles and low hydraulic roughness.The mixed tidal-fluvial regime is characterized by a weaker tidal influence, a shallower and more variable water depth, lower maximum absolute shear stresses, asymmetric dunes, higher leeside and slip face angles, and a rougher hydraulically regime.The increase in depth is the main reason that hydraulic roughness is lower in the tidal regime (Equation 13), since the sources of roughness in the Main Channel, sediment composition and dune height, are relatively constant.Contrasting flow conditions in the two regimes are not reflected in dune height or length.In other systems, dune height is sometimes found to decrease in tidally influenced regions (Prokocki et al., 2022).Rapid local deposition of the sediment in the deltaic part of the Fraser might result in tidal dunes that are larger than expected (Villard & Church, 2005).The change in flow regime is reflected in the leeside angle, slip face angle, dune symmetry and dune crest shape.In particular slip face angles are significantly larger in the fluvial-tidal regime, on average 13° compared to 7° in the tidal regime.Dunes are on average asymmetric upstream of the bifurcation at RK 40 (Figure 7), and symmetric downstream of RK 40.This agrees with the findings of Kostaschuk and Villard (1996), who relate the symmetric dunes to high sediment transport rates due to the tides, potentially flattening the leeside angles.Indeed, high maximum shear stresses (Figure 11b) are observed in the tidal regime, although tidally averaged shear stresses remain relatively constant (Figure 4b).
The reversal of the current switches the leeside and stoss side every tidal cycle, steepening both sides in a similar manner (Lefebvre & Winter 2016).This could be one of the reasons for the large observed variability in angles, since the MBES data is simply a snapshot of the riverbed.Bidirectional currents cause crest orientation to be time-dependent (Hendershot et al., 2016).Both the duration (t rev ) and the strength of the flow reversal (τ or Q) determine the dune shape.During low river discharge conditions, the maximum upstream-oriented discharge at RK 22 varies between 4,000 and 6,000 m 3 s −1 , depending on the spring-neap tide cycle.Only 30 km further upstream this has decreased by 66%-75%, although the reversal time has only dropped by 9%.(fluvial-)tidal dunes, and fluvial dunes.The former were restricted to the most downstream reach (RK < 30 km), and were upstream oriented, predominately low-angled (based on maximum LSA), 2D dunes.Fluvial dunes were downstream oriented, and were higher and longer than the tidal dunes, with higher slip face angles.The division of the regimes in the Columbia is clearer than in the Fraser.One of the reasons for this could be that the division in the Columbia coincides with a change in grain size.The relatively minor downstream fining in the thalweg of the Fraser river (Figure 2e), which this is not uncommon in other large sand-bedded rivers (Frings, 2008), might be one of the contributors to the absence of a decrease in dune height in downstream direction as observed in the Lower Columbia.Lefebvre et al. (2021) also found an increase in dune symmetry in the downstream direction in the Weser Estuary, Germany, but they did not distinguish between two different regimes.However, their data shows that around 60 km from the river mouth, upstream of the estuarine turbidity maximum, the leeside angle of dunes decreases, and dunes become more symmetric.The gradual transition is almost twice as far upstream as in the Fraser River, potentially attributable to the tidal effect in the Weser extending much further upstream than in the Fraser.In this study, and in those of Prokocki et al. (2022) and Lefebvre et al. (2021), the transition in dune morphology coincides with an increase in channel cross-sectional area, either by widening, deepening or both.The deeper regimes are more tidally dominated, and the constriction upstream of the cross-sectional area leads to a rapid dissipation of tidal energy, that is reflected in the dune leeside angle and symmetry.Comparison between different river systems and disentangling the drivers behind the different patterns in dune geometry is worth future study.

How Can Dune Variability in the Fluvial-To-Tidal Transition Zone During Low River Discharge Be Estimated and Explained?
Tidally averaged bed shear stresses from the hydrodynamic model can be used to reliably estimate reach-averaged dune height using the predictor of Van Rijn (1984).Furthermore, the shear stress distribution estimated by the hydrodynamic model with constant roughness can be used to estimate local dune patterns (Figures 8b-8g), thereby capturing the cross-sectional variability in dune heights as observed in Figure 6.Shear stress, which varies over the cross-section, and is one of the input parameters of the dune predictor of Van Rijn (1984), largely explains the observed patterns.For example, in zone 1, dune height decreases downstream, because river width increases and flow velocity decreases, resulting in lower shear stresses (Figure S10 in Supporting Information S1).In zone 2, dunes are highest on the south side of the channel, where the river is deepest and the flow velocity and shear stresses are highest.Finally, in zone 3, centrifugal acceleration generates higher flow velocity and larger dunes on the outside of the bend, whereas upstream the dunes are the largest on the inner bend because, flow is accelerated by the momentum inherited from the bend upstream (Jackson, 1975).
Van Rijn (1984) and other dune height predictors tested did not yield accurate estimates of absolute magnitude of local dune height using tidally averaged bed shear stresses from the hydrodynamic model.Therefore, we did not use the predicted values in the roughness predictor, but used averaged values instead.The dune height predictors do reproduce the overall patterns of dune height, suggesting that the right processes are captured by the predictors.The Van Rijn (1984) calibration factor may be adjusted to better represent the studied field conditions.The poor estimation of absolute values is likely related to a number of factors that are not included in the predictors, including self-organization dunes in a shear stress field (Bradley & Venditti, 2019) (such as merging and splitting (Gabel, 1993;Hendershot et al., 2018), crest line deformation (Venditti et al., 2005)), local sediment dynamics not captured by low resolution bed sediment sampling such as local scour (Leclair, 2002), discharge fluctuations and associated hysteresis (Bradley & Venditti, 2021;Julien et al., 2002) and the potential presence of remnant dunes from earlier high-flow conditions.
The influence of local factors can be seen in our three focus zones.The larger dunes observed in zone 2 may be related to the local sediment supply being higher here, so dunes develop to a maximum equilibrium size compared to zones 1 and 3.The dunes in zone 2 become longer in the downstream direction until they disappear, even though flow velocity and grain size do not change significantly.The disappearance of dunes in this area could be because the surface of the bed consists of a thin layer of medium sand overlying a deposit of Pleistocene or early Holocene sediment, such as cohesive clay (Clague et al., 1983) (see Figure S1 in Supporting Information S1), that is not conducive to dune formation.Similarly, in zones 1 and 3, dunes do not develop where the outer bank cuts into a clay layer (Figure S9 in Supporting Information S1).Additionally, the dunes could be reworked remnants from the previous freshet (Bradley & Venditti, 2021), and their geometry could be related to the much stronger and predominantly downstream currents associated with high river discharge.However, Kostaschuk et al. (1989) found that dunes near Steveston (RK 10) adjusted to the post-freshet decline in discharge over a period of weeks, supporting our contention that the dunes observed herein (more than 6 months after the last freshet) are at least in quasi-equilibrium with low-flow conditions.Bradley and Venditti (2021) interpreted low-amplitude bed undulations at RK ∼35 as relics from higher flow conditions with smaller dunes superimposed, the latter formed by the low-flow conditions.Kostaschuk et al. (1989) interpreted similar features as "washed-out" dunes that represented a transition from large, freshet bedforms to small dunes adjusted to low river discharge.In this study, we detrended the bed level prior to measuring dune geometry, thereby ensuring that the dunes that we analyzed were representative of low flows.
The poor estimation of local dune geometry in the FTTZ and the observed variability in dune morphology has practical implications for scientists and engineers.First, fairway depth cannot be maintained solely on the basis of an average dune height, because height varies unpredictably over the river bed.Second, models based on reach-averaged dune geometry may result in inaccurate estimates of form roughness and water levels and local values should be used instead.Third, using estimated dune height to calculate hydraulic roughness would result in an inaccurate roughness prediction in the FTTZ.Finally, because the variability in dune height across the channel is twice that of dune height variation along the channel, the grid cell size in hydraulic models should be twice as large in the longitudinal direction than in the cross-river direction.

To What Extent Does Dune Geometry and Variability Affect Reach-Scale Hydraulic Roughness?
Hydraulic roughness is traditionally predicted using dune height and length (Bartholdy et al., 2010;Lefebvre & Winter 2016;Soulsby, 1997;Van Rijn, 1984).However, recent research shows that the leeside angle of dunes might be important for roughness prediction (Lefebvre & Winter 2016).Characteristics of the leeside angle determine the strength of the flow separation zone (Lefebvre et al., 2014) in turn affects form roughness (Lefebvre et al., 2013).Large rivers are covered by low-angled dunes (LADs) with slip face angles (SFAs) < 30° (Cisneros et al., 2020;Kostaschuk & Venditti, 2019) that generate less flow separation than high-angled dunes that have steeper SFAs (Kwoll et al., 2016).However, weak or intermittent flow separation, with mean leeside angles (LSAs) of only 10° still generates flow resistance (Kwoll et al., 2016).Lefebvre and Cisneros (2023) show that not only the leeside angle itself, but also the shape of the leeside impacts flow properties and turbulence.Based on numerical experiments, they found that LADs with a mean LSA of <10° and a SFA of <20° are not capable of permanent flow separation.LADs do generate turbulence (Kostaschuk & Villard, 1996;McLean & Smith, 1979), however, because the decelerated downstream flow generates a shear layer that causes eddy generation (Best & Kostaschuk, 2002;Kostaschuk & Villard, 1999), sand resuspension and roughness.
In this study, the transition from a fluvial-tidal to a tidal regime and the corresponding change in dune LSAs and SFAs are not reflected in the reach-scale hydraulic roughness needed to attenuate the tidal motion in the model.Implementing a roughness change at the depth break at RK 40 could parameterize the change in dune leeside angle at the regime transition.Models with a higher roughness downstream than upstream performed slightly better than models with the highest roughness upstream (Text S3 in Supporting Information S1).This is contrary to the expectations based on the leeside angle observations, and suggests a different source of roughness in the tidal regime (see below).Additionally, the dune roughness predictor (Equation 2), based on dune height and dune length, yields very similar values to the calibrated model roughness (RMSE f = 0.0053).This supports the application of the dune roughness predictor in a tidal environment, and also indicates that dune leeside angle might not be important in determining reach-scale model roughness.Finally, hydraulic model performance is not improved by using local dune geometry.This in turn suggests that variable dune roughness might not be needed to simulate large-scale (tidal) flow.Similar conclusions were drawn from a fluvial system where dune roughness calculated from dune geometry explained only 31% of the variance of the roughness inferred from the water surface slope (de Lange et al., 2021) and the remaining variance could not be explained by leeside angle statistics.
The limited impact of variation in local dune height, length and leeside angle on the hydraulic model could be due to several factors.First, 3-dimensional dune fields such as those in our study area, generate less roughness than 2-dimensional dune fields (Venditti, 2007), which could explain the lack of model improvement when implementing dune-related roughness.Second, a complex leeside shape might have an effect on flow separation (Lefebvre & Cisneros, 2023) and form roughness.So, even though the SFAs are large enough to generate flow separation, the shape of the leeside might prevent it.Third, we evaluated the hydraulic model by assessing tidal attenuation and water level fluctuations and found minimal impact from local dune geometry.This notwithstanding, incorporating dune roughness could be important for estimation of residual sediment transport (Brakenhoff et al., 2020;Herrling et al., 2021).Local values of shear stresses may be required for morphodynamic modeling, even if they are not needed for modeling tidal propagation in a hydrodynamic model.
We evaluated the model on the reach-scale where other components of roughness dominate (see below).These components are less relevant on the local scales, where dunes might be the main source of roughness.In addition, the estimation of hydraulic roughness generated by dunes deviates locally from the constant model roughness.
As a result, in the mixed tidal-fluvial regime, the dune-induced roughness is larger than needed for attenuation based on the calibrated roughness.Davies and Robins (2017) found that the overall effective roughness of the bed is about half of the maximum local dune-induced roughness expressed in k s .Halving the k s value in Equation 13results in a comparable dune roughness and calibrated roughness in the mixed fluvial-tidal regime (RMSE 0.0034 for RK > 40; Figure S4 in Supporting Information S1), but not in the tidally dominated regime where dune roughness remains lower than calibrated roughness.This could be a result of lower LSAs in the tidal regime.
Including the LSA in roughness estimation using the equation developed by Lefebvre and Winter (2016) results in unrealistically low values of roughness.In general, the available evidence indicates that the LSA does not have a significant impact on reach-scale model roughness.
In our research area there are several reach-scale sources of roughness.First, large-scale river geometry, which is included in the hydraulic model by the bathymetry.We observed an out-of-phase relation between hydraulic roughness and the smoothed gradient of the bed level in the tidally dominated regime (Figure 10).A similar relation was observed in the Rhine and Waal rivers in the Netherlands by de Lange et al. ( 2021), who hypothesized that multi-kilometer depth oscillations induce flow divergence associated with depth increase, which in turn causes energy loss.This is reflected in an elevated hydraulic roughness.Such an out-of-phase relation is again seen in the tide-dominated part of the lower Fraser, but not in the mixed tidal-fluvial regime.There, depth increases coincide with decreases in width, keeping the cross-sectional area relatively constant.As a result, changes in depth do not result in flow divergence or convergence and the out-of-phase relation does not develop.Second, intertidal areas affect reach-scale roughness.The calibrated friction in our model is an indication of the friction required to attenuate the tide.The model calibrated with a uniform Manning's roughness (n man = 0.026 s m −1/3 ) performs reasonably well in modeling tidal propagation, but regions with a significant decrease in tidal energy (between RK 9-18.5 and 35-42) are not well-captured by the model (Figure 3).These regions include intertidal areas (Figure S11 in Supporting Information S1), which flood during high tide and are not properly represented in the model, limiting tidal attenuation.Also, engineering works, such as the tunnel at RK 18 and the bridge at RK 36, could be an extra source of roughness.
It is worth noting that the results from this study may not be directly applicable to higher flows in the Fraser or other FTTZ-systems.During high discharge conditions in the Fraser (freshet), dunes generally increase in size compared to low flows (Bradley & Venditti, 2021), and could be a larger source of roughness.However, slip face angle does not seem to increase significantly during freshet conditions (Bradley & Venditti, 2021), meaning that flow separation and the resulting elevated roughness will not occur more frequently than during low discharge.Additionally, during the freshet, more large-scale roughness elements will be included in the river corridor, including submerged intertidal areas and bars.We consider the effects of large-scale geometry changes on hydraulic roughness to be a key subject that has remained severely understudied to date.

Conclusions
About 40 km from the mouth of the Fraser River, a local depth decrease separates a fluvial-tidal regime at the landward side from a tide-dominated regime crossing the Fraser Delta.Bathymetric data and a hydraulic model of the lowermost 85 km of the river were used to explore the spatial variability and controls of dune morphology in this fluvial-to-tidal transition zone (FTTZ).From our investigations, we conclude that: • There are no significant spatial trends in dune height or length, even though the river deepens at RK 40.Local variability in dune height and length dominates regional trends, and variability in dune height and length in the cross-sectional direction is twice as large as in the longitudinal direction.• In the tidal river, dune height predictors provide a good first approximation of regional dune height and local spatial patterns, but local shear stress predictions need to be improved to enable robust estimates of local dune heights.Using shear stresses from the hydraulic model calibrated with a constant Manning's n roughness of 0.026 s m −1/3 , the dune height predictor of Van Rijn (1984) reproduces local patterns of dune heights, using tidally averaged values of bed shear stress.Other predictors of dune height perform worse.• Mean leeside and stoss side angles of dunes are lower in the tidal regime compared to the fluvial-tidal regime, and dunes become symmetric under tidal influence.These changes in dune morphology do not affect reach-scale hydraulic roughness, because the calibrated model roughness is similar to the dune-derived roughness based on dune height and dune length.As a result, hydraulic model performance using a calibrated, constant, roughness is not improved by implementing dune-derived bed roughness.• In the Fraser River, large-scale variations river geometry are more important than dune morphology in controlling variations in reach-scale roughness, a finding which is not unique to the Fraser River (de Lange et al., 2021).In the Fraser river, oscillations in river depth elevate hydraulic roughness in the tidal region, but this does not occur in the fluvial-tidal regime because changes in depth are compensated for by changes in width, keeping the cross-sectional area of the channel relatively constant.Intertidal areas in the Fraser are likely a significant source of flow resistance.Flow resistance by depth variation and intertidal areas needs further study.
(a)   How are bedform characteristics impacted by the sudden change in tidal flow strength?(b) How can dune variability in the fluvial-to-tidal transition zone be explained and estimated?(c) To what extent does dune geometry and variability affect reach-scale hydraulic roughness, and which other factors can play a role in determining this bed roughness?Bathymetric field data from base flow conditions were used, allowing us to focus on the impact of the tides, which penetrate further upstream during base flow.A 2D hydrodynamic model is created to assess hydraulic roughness, and to explore the impacts of spatial variation in river and tidal flow.The field data is analyzed to assess bedform characteristics and its relation to hydraulics and river geometry, and model output is used to estimate dune geometry and assess its impact on hydraulic roughness.

Figure 1 .
Figure 1.Study area of the Fraser River in British Columbia, Canada (a).(b) The Fraser River from river kilometer 10 (Steveston Harbor) to 85 (Mission).Green shaded area indicates the model domain.Gray markers indicate gauging stations.(c-e) Measured multibeam bathymetry in the three focus zones examined in this study, shown in hill shade, (f-h) examples of hill shade depiction of the measured dune fields.(i-k) example profiles of the dune fields along the corresponding black lines in subfigures (f-h).

Figure 2 .
Figure 2. Grain size distributions of bed sediment in the Fraser River, shown in distance from the river mouth at RK 0 (chainage).(a-c) grain size distribution along the north bank, thalweg and south bank.(d) cumulative distribution at the north bank, thalweg and south bank.(e) median grain size (D 50 ) in and outside of the thalweg.Markers differentiate between samples taken in the thalweg (solid marker) or outside along the river banks (indicated with "x").The data is from a data compilation byVenditti and Church (2014), which includes reanalyzed observations fromMcLaren and Ren (1995), and recent observations byMurphy (2023).

Figure 3 .
Figure 3. (a) Calibration of the model with uniform roughness, shown along the river (chainage, where 0 indicates the river mouth).The observed tidal amplitude of the tidal constituents M2 (black bars), K1 (dark gray bars), and M4 (light gray bars), and the corresponding modeled tidal amplitudes are indicated.(b) Discharge at Mission.Highlighted part of the discharge curve indicates the timeframe of MBES data collection.(c) modeled water surface slope over time, simulated with the model with n man = 0.026 s m −1/3 .(d) modeled propagation of the tidal wave per station, simulated with the model with n man = 0.026 s m −1/3 .

Figure 5 .
Figure 5. Morphological characteristics of the Lower Fraser River along chainage (river kilometer).(a) smoothed channel width (W) and smoothed cross-sectional area (A), (b) bed level z in meters above sea level, (c) channel shape, expressed as depth excess (D e ), transverse bed slope (ξ) and curvature (r).

Figure 6 .
Figure 6.Dune geometry along chainage (river kilometer).(a-c) Dune height (Δ; black) and dune length (λ; blue) throughout the research area.Human-made structures, dredging marks, confluences, bifurcations and bars, focus areas, and zones with no data are indicated (see legend).(d, e) Standard deviation (σ) within the mean multibeam echosounder coverage width (230 m) of dune height and dune length over the cross-section, north bank, thalweg and south bank.In each bar, the central mark indicates the median, and the bottom and top edges of the box indicate the 25th and 75th percentiles, respectively.The whiskers extend to the most extreme data points, and outliers are not shown.

Figure 7 .
Figure 7. Dune shape in the study area along chainage (river kilometer).(a) leeside angle (LSA) and stoss side angle (SSA), (b) dune slip face angle (SFA), (c) dune asymmetry, expressed as the ratio between the length of the (seaward) leeside and the total bedform length.A value of 0.5 indicates symmetric dunes, values of asymmetry smaller than 0.4 are defined as flood-asymmetric, while values larger than 0.6 are ebb-asymmetric.Confidence intervals of linear regressions are shown.Subfigure (d, e) show dune leeside angle against reversal time (t rev ) and against maximum shear stress (τ max ), respectively.

Figure 8 .
Figure 8. Residual dune height (measured dune height Δ meas minus estimated dune height Δ est ) along chainage (river kilometer), to assess dune height predictor performance of the predictor of (a) Van Rijn (1984), (b) Yalin (1964), (c)Karim (1995), and (d)Venditti and Bradley (2022).The measured data is based on the average of three longitudinal transects, and includes the minimum and maximum values in a gray shading.The modeled data is based on the model with a constant roughness of n man = 0.026 s m −1/3 .In subfigure (a) estimated dune height with dune-adjusted roughness (varying between 0.024 and 0.028 s m −1/3 ) is displayed in green.

Figure 10 .
Figure 10.Hydraulic roughness along the river (chainage).(a) smoothed roughness (f) (and original roughness in gray) calculated from dune geometry (Equation 13)(blue) and the gradient of the smoothed bed level (∂z/∂s; black).(b) roughness expressed in Manning's roughness coefficient n man .Model roughness with a constant n man of 0.026 s m −1/3 (black), roughness estimated from dune geometry (Equations 13 and 17) (blue), and dune-adjusted model roughness (green).