Spatio-temporally Varying Manning Roughness in Rivers and Streams: A calibration approach using in-situ water level and UAS altimetry

11 Hydraulic roughness (expressed in terms of e.g. Manning's roughness coefficient) is an important input to hydraulic 12 and hydrodynamic simulation models. One way to estimate roughness parameters is by hydraulic inversion, using 13 observed water surface elevation (WSE) collected from gauging stations, satellite platforms or UAS (Unmanned 14 Aerial System) -based altimeters. Specifically, UAS altimetry provides close to instantaneous observations of 15 longitudinal profiles and seasonal variations of WSE for various river types, which are useful for calibrating 16 roughness parameters. However, it is computationally expensive to run high-resolution hydrodynamic models for 17 long simulation periods (e.g. multiple years), and thus global optimization of spatially and temporally distributed 18 parameter sets for such models, e.g., spatio-temporally varying river roughness, is still challenging. 19 This study presented an efficient calibration approach for hydraulic models, using a simplified steady-state 20 hydraulic solver, UAS altimetry datasets, and in-situ observations. The calibration approach minimized the 21 weighted sum of a misfit term, spatial smoothness penalty, and a sinusoidal a priori temporal variation constraint. 22 The approach was first demonstrated for several synthetic calibration experiments and the results indicated that 23 the global search algorithm accurately recovered the Manning–Strickler coefficients M for short river reaches 24 in different seasons, and M varied significantly in time (due to the seasonal growth cycle of the aquatic 25 vegetation) and space (due to, e.g. spatially variable vegetation density). Subsequently, the calibration approach 26 was demonstrated for a real WSE dataset collected at a Danish test site, i.e., Vejle Å. Results indicated that 27 spatio-temporal variation in M was required to accurately fit in-situ and UAS altimetry WSE observations. This 28 study illustrated how UAS altimetry and hydraulic modeling can be combined to achieve improved 29 understanding and better parameterization of small and medium-sized rivers, where conveyance is controlled 30 by vegetation growth and other spatio-temporally variable factors. 31

The approach was first demonstrated for several synthetic calibration experiments and the results indicated that 23 the global search algorithm accurately recovered the Manning-Strickler coefficients M for short river reaches 24 in different seasons, and M varied significantly in time (due to the seasonal growth cycle of the aquatic 25 vegetation) and space (due to, e.g. spatially variable vegetation density). Subsequently, the calibration approach 26 was demonstrated for a real WSE dataset collected at a Danish test site, i.e., Vejle Å. Results indicated that 27 spatio-temporal variation in M was required to accurately fit in-situ and UAS altimetry WSE observations. This 28 study illustrated how UAS altimetry and hydraulic modeling can be combined to achieve improved 29 understanding and better parameterization of small and medium-sized rivers, where conveyance is controlled 30 by vegetation growth and other spatio-temporally variable factors. geometry (which may only be available for a few points along the river), is a central task in the development of 48 real-world hydraulic and hydrodynamic models. River roughness is an effective parameter representing friction 49 effects in the shallow water equations. It is a critical controlling parameter for conveyance estimation and thus 50 water levels in rivers and streams. The parameter is influenced by many factors, such as the type of the bed and 51 bank materials, aquatic vegetation, surface irregularity, shape and size of the channel cross-section, the 52 meandering character of the river channel, and state of flow motion (Cowan, 1956). Some studies have 53 developed a predictive model for mean flow in irregular natural rivers, and the results indicated that the effective 56 resistance was strongly influenced by river cross-sectional nonuniformity, and the authors pointed out that 57 sampling density for geometric parameters should depend on the degree of stream irregularity. Tuozzolo et al. 58 (2019) further analyzed the impact of reach averaging Manning's equation and showed that roughness varied 59 significantly even in a 6.5 km river stretch. 60 Meanwhile, seasonal variations in discharge and aquatic vegetation (type, height, and density) also 61 profoundly impact the flow resistance, especially for many lowland vegetated river channels (Marjoribanks et 62 al., 2014). The seasonal growth of aquatic macrophytes can significantly increase bed resistance, which leads 63 to a decrease of river channel conveyance and consequently may increase flood risk. Aquatic vegetation in rivers the Gauckler-Strickler coefficient Ks (Ks = 1/n) increased significantly after vegetation was cut. However, 66 simultaneously calibrating roughness parameters in both spatial and temporal scales is seldom reported. 67 The calibration of spatiotemporally varying river roughness requires dense sampling of the water 68 surface elevation (WSE) in both time and space. Gauging stations and satellite altimetry are two widely used 69 methods to retrieve WSE observations. Gauged WSE is only available at discretely distributed points in space, calibrating parameters of hydrodynamic models due to the high computational cost of the forward simulation 88 models. However, such algorithms cannot guarantee a global optimum and may perform poorly in highly 89 parameterized and non-linear inverse problems. Global optimization methods are more suitable and have a high 90 probability of successfully finding the optimal parameter set (Duan et al., 1992(Duan et al., , 1993. Still, global search 91 algorithms require more runs of the forward simulation model (Kittel et al., 2021). 92 This study proposed the estimation of optimal spatiotemporally varying Manning-Strickler coefficients 93 using station data and UAS altimetry. Considering station data characterized by high-temporal resolution (daily 94 or sub-daily) but coarse and discrete distribution in space, and UAS altimetry with high-spatial resolution and 95 accuracy (centimeter level) but sporadic coverage in time, the present study explored the value of 96 simultaneously using both in-situ and UAS observations of WSE for Manning-Strickler coefficients calibration. 97 The methodology was first evaluated using different synthetic calibration experiments, considering data roughness parameters for the stream was calibrated and validated using real-world WSE data collected from 100 We assumed that the effective roughness in each small river stretch x was constant, and the roughness 110 parameters were also stable over a short time interval t, e.g., a month in this study. Thus, we had a 111 spatiotemporally varying-parameter M(x, t), the symbol M is the Manning-Strickler coefficient; M = 1/n, n 112 being Manning's roughness coefficient. Considering that a large number of parameters needs to be calibrated, 113 the objective function for the optimization algorithm should be regularized. Here, the objective function 114 comprised a !"#$"% term and two '()*+,'"-,%"./ terms. The !"#$"% term was set to measure the weighted 115 squared difference between the simulated WSE and the observed truth. The '()*+,'"-,%"./ terms were used 116 to constrain the variations of the parameters. The objective function (∅) was presented as: 117 ∅ = 2 • !"#$"% + (1 − 2) • '()*+,'"-,%"./ (1) 118 Where 2 was the weight between !"#$"% , and '()*+,'"-,%"./. 119 The misfit is the difference between the model simulations (9:; <=> ) and the observations (9:; ?@< ) 120 at all the time steps (A B ) and chainage points (A C ). The observations included in-situ WSE and UAS altimetry, Where O is the weight between smoothness and a priori term. 126 The smoothness term was based on the fact that the river channel conditions (e.g., size and type of the continuously in space and the roughness parameters along the river channel was expected to change smoothly. 129 The smoothness term was expressed using a standard first-order roughening term, i.e. 130 Where the D <>??BU is a typical variation between neighbors in space. 132 Considering that the inter-annual and seasonal variations of the discharge and vegetation growth periods, 133 we assumed that the temporal variations of river roughness have a sinusoidal regular pattern. For example, the Where D YZV=?V= is a typical deviation from the a priori sinusoidal seasonal model.  Thus, an efficient steady-state hydraulic solver was used. 159 The simplified hydraulic solver is based on the Saint-Venant equations, assuming that the water flow 160 is in the condition of steady-state, i.e., all quantities are invariable in time. The development of the solver is 161 described fully in Kittel et al. (2021). Only the main formulas of the solver therefore are presented here. 162 The momentum balance equation under steady-state is presented as: 163 Where A is flow cross-sectional area, Q is volumetric discharge, x is the chainage, g is the gravitational 165 constant, h is flow depth, : [ is the bed slope, : j is the friction slope. The general form of the equation was 166 organized by taking the partial derivative of the area relative to the chainage and width, and isolating the change 167 in depth over the chainage: 168 Where : [ is calculated using river bed elevation (z) and chainage as ( − g-gG ⁄ ), : j is the river channel R is the hydraulic radius. Equation (9) was solved explicitly. We used a grid spacing of ΔG =20 m, and 173 the hydraulic parameters, i.e., area, depth, bed elevation, were interpolated for each calculation grid point from 174 neighboring cross sections. alternating q-and h-grid points, i.e. points where the discharge, q, and water level h, respectively. The 182 computational grid is generated automatically based on the user requirements, but the lower value of the grid 183 spacing needs more calculations.
Vejle Å (Vejle River) is a natural stream located in Vejle Municipality, Denmark, which is a typical 186 lowland river with mild slopes and meandering character. It originates from Engelsholm Sø and empties into 187 Vejle Fjord, with an approximate length of 36.7 km (Fig. 1). The stream flows through woods, farmland, bush, 188 grassland, and urban regions, with the riverbank densely covered with shrubs, trees, and weeds. The dense 189 vegetation in the Vejle Å increases friction in the river channel and raises the water level. It is necessary to cut 190 the submergedvegetation in the river during the flooding season to improve river channel conveyance.  The annual average discharge in Refsgårdslund and Haraldskaer were 3.03 m 3 /s and 4.19 m 3 /s, respectively. The maximum daily discharge during the period is 10.68 m 3 /s and 16.13 m 3 /s. There are no influential tributaries in 207 the chainage interval between Refsgårdslund and Haraldskaer, and, therefore, we assumed that the lateral inflow 208 is uniformly distributed between stations Refsgårdslund and Haraldskaer. 209

Water surface elevation 210
Water surface elevation data was collected from five hydrological stations, i.e., Tørskind (14.77 km), 211 Ravning (18.92 km), Vingsted (21.29 km), Rosborg (33.56 km), and Vejle Havn (36.57 km), and the data are 212 displayed in Fig. A2. The time resolution of the WSE data is 15 minutes. The river surface has an average slope 213 of 5.8 ‰ between upstream Refsgårdslund and the outlet to the sea. WSE of Vejle Havn was used as the 214 downstream boundary of the hydraulic and hydrodynamic models while the other stations were used for 215 calibration. 216

Cross sections 217
The cross-section information was obtained from the national information system run by WSP Denmark, 218 and 130 cross-sections were included from Refsgårdslund to Vejle harbor. The raw cross-section data includes 219 the distance from the left levee and the corresponding depth. The raw data were then put into MIKE Zero 220 software for further processing. The processed cross-sections contain hydraulic parameters including width, 221 hydraulic radius, bed elevation, submerged cross-sectional area, and the depth at a specific WSE. The basic 222 information of all bank-full cross-sections is shown in appendix (Fig. A3). The river channel is narrow and deep 223 in the upstream area, wide and shallow in the downstream area. The bank-full submerged area does not increase 224 from upstream to downstream, which is a likely cause of increased flood risk, e.g., the cross-sectional area near 225 Haraldskaer located at chainage of 28.13 km is reduced and this area suffered from frequent floods. The river 226 The UAS-borne WSE measurement system used in this study was developed by Bandini et al. (2020). 230 This system can measure accurate and distributed high-resolution WSE for small-size and vegetation-covered 231 streams. The system uses a lightweight hexacopter drone equipped with a dual-frequency differential Global 232 Navigation Satellite System (GNSS) and a 77 GHz radar chip with full waveform analysis to measure WSE 233 (Bandini et al., 2020). 234 Several field campaigns were implemented from 2018 to 2020 by DTU Environment and Drone 235 Systems (https://dronesystems.dk). Two long-distance UAS surveys were carried out in February and June 2020, 236 which covered a river reach with a total length over 20 km (Fig. 2). The weather conditions and power supply 237 infrastructure are the main constraining factors for the drone survey, and a one-kilometer-long river stretch can 238 be measured in approximately 15 minutes. UAS altimetry data processing included post-processed kinematic processed data was further filtered by a river mask to guarantee that the radar observations were captured above 241 river water surface. The water mask was generated by buffering the river center line by 5 meters in each side. 242 Detailed information of the UAS radar altimetry post-processing can be found in Bandini et al. (2020). 243 UAS altimetry shows great potential in delivering high-resolution WSE observations, as illustrated in 244 Manning-Strickler coefficients M is created with seasonal sinusoidal variation (Fig. 3). 261   is 0.46 m 1/3 /s for Scenario 2 with sufficient WSE. However, the WSE was over fitted by Scenario 2 (Table 1), 324 which may indicate that a larger weight for the regularization terms may be appropriate. The RMSE of calibrated 325 M is 2.53 m 1/3 /s for Scenario 3, which is inferior to the calibration results of Scenario 2 due to the limited amount 326 of available WSE. The values of M lower than 20 m 1/3 /s were significantly overestimated because the available 327 WSE observations are insufficient for adequately constraining the parameter. 328 The evaluation results of the simulated WSE against synthetic truth in temporal and spatial scale showed 329 that the accuracy of the simulated WSE is high in both time and space under Scenario 2 and Scenario 3. In 330 Scenario 2, even though the RMSE shows some spatial (red columns in Fig 5e) and temporal (gray columns in 331 the Fig 5f) difference, the overall RMSE is low. The RMSE of the simulated WSE is high in Scenario 3, which 332 indicates that the available WSE is insufficient to constrain the variations of M. The RMSE is high in low-flow 333 seasons, e.g., July to December, but the RMSE is low in high-flow seasons, i.e., January to May. Along the chainage, the simulate WSE is less accurate for chainage 20 -30 km, because M was overestimated compared 335 to the synthetic truth in this chainage interval.

Real-world calibration for specific dates 344
The synthetic studies demonstrated the ability of the calibration approach to fit spatiotemporally 345  The calibration results of the M and the comparison between simulated and observed WSE are displayed 372 in Figure 7. The spatio-temporal variation of M is distinctive, high in the cold season (December to April) and 373 low in the warm season, relatively high upstream of the river and low downstream. In the upstream Vejle Å (12 374 -24 km), the bed elevation shows undulating terrain, and the bankfull width and area are changeable in this 375 area compared to the downstream (Fig.A3). Thus, the complex bathymetry coupled with seasonal growth of 376 vegetation, resulted in significant changes of M in different seasons in this area. The calibrated M was more 377 inconsistent in cold season due to the high-flow but relatively stable in warm-season with low-flow. UAS 378 altimetry and in-situ observations of WSE were used for the steady-state solver calibration, and the RMSE of 379 the simulated WSE is 7.67 cm (table 2). While a significant bias between simulated WSE and UAS altimetry is 380 visible, as shown in Fig. 7b. Because a weight between drone data and stations data was used to balance their 381 corresponding misfit in the objective function. The results indicate that the weight of in-situ data is too large.  The comparison results between the simulated WSE against in-situ observations of four hydrological 390 stations are displayed in Figure 8. Generally, the simulated WSE agreed well with in-situ observations in the 391 high-flow period, and the simulated WSE is clearly less accurate in the low-flow seasons, which is similar to 392 the synthetic study. The overestimation of WSE in low-flow seasons indicated an overestimated of M. 393 Considering the UAS altimetry data we used for calibration in different seasons, the WSE on 02/24/2020 used 394 sufficiently for constraining the M. 395 Validation results indicated that the calibrated M were not representative for the validation period, 396 especially at Vingsted station (Fig. 8). Spatial patterns and dynamics of vegetation growth may change from 397 year to year in which case a longer calibration period is required to obtain more robust parameter estimates. 398  In this study, we used a radar altimetry payload on a UAS to measure the WSE of Vejle Å. The altimetry 418 data was smoothed with an interval of 20 m, equivalent to the grid resolution of the steady-state hydraulic solver, 419 ensuring that at least one WSE observation fell on each calculation node. The results indicated that the 420 parameters were well constrained by dense WSE datasets and the uncertainty had significant effects on the 421 calibration results, as illustrated by the synthetic experiments. Potentially, we could further increase the grid 422 resolution to 1 m or even less, considering the spatial sampling resolution of WSE by UAS. However, 423 calculation efficiency will be decreasing as the number of grid points increased.
The advantages of using UAS altimetry data along the chainage were also demonstrated by comparison 425 to the calibration experiments in which temporally continuous in-situ observations and the MIKE Hydro River 426 model were used. The calibrated results improved slightly after including UAS altimetry (Fig. B1). Additionally, 427 UAS altimetry data can be used to calibrate highly resolved Manning-Strickler coefficients M for a particular 428 period of interest, for example a flooding period, which was not possible when using in-situ station data only. 429 Additionally, the coverage of the UAS altimetry is essential when spatially variable roughness 430 parameters are calibrated. As shown by Scenario 1, the calibration results of M in 06/02/2020 outperformed 431 02/24/2020. Available UAS altimetry were more continuous in 02/24/2020; even though the spatial variations 432 of M were constrained, the fitted roughness coefficient was far away from the truth in the chainage where the 433 UAS altimetry was missing. In contrast, the available UAS altimetry was more scattered on 06/02/2020. 434 Although UAS altimetry was missing for short intermittent reaches, M can be estimated in those reaches thanks 435 to the smoothness regularization term. However, M was not well-constrained for longer reaches without 436 available UAS altimetry data example can be seen in the results for 06/02/2020. We can thus conclude that for The calibration approach for spatial-temporally varying M can also be used for the rivers with complex 453 topography and changeable hydraulic conditions. For instance, M may vary significantly at high stage for 454 braided and terraced river. Meanwhile, river resistance will be significantly increased by river ice-jams that 455 cause severe floods in some middle and high latitude rivers, e.g., Yellow River. Thus, a varying M might 456 increase the performance of hydrodynamic simulations in such situations.
The steady-state solver is efficient and effective for calibrating the M as shown in the present study. It 459 took millisecond for the Steady-state solver to simulate WSE in Vejle Å with a resolution of 20 m. Compared 460 to the hydrodynamic MIKE HYDRO River, the steady-state solver is computationally efficient. To go a step 461 further, we can use global optimization to fit M by using this solver. The global optimization is independent of 462 initial value and capable of avoiding local optima. To calibrate the spatiotemporally distributed M, the number 463 of parameters increased significantly. The high resolution of the parameters in space and time increases the 464 computational requirements for streams or large rivers (Kittel et al., 2021). Moreover, calibration results will 465 be significantly affected by the chosen regularization strategy. 466 The present calibration algorithm for spatio-temporally distributed M has several drawbacks. In the 467 hydraulic model, we used the cross-sectional information to delineate the river channel, and the cross-sections 468 were linearly interpolated for each calculation node from the limited in-situ data. Thus, the accuracy, 469 representation and density of the surveyed bathymetry is critical for the calculation. The 1-D steady-state 470 hydraulic solver is based on the Saint-Venant equation, and it is insufficient to simulate overbank water flow. 471 For such situations, a 1D-2D hydraulic model would be required. That is a reason for the relatively high error 472 of the simulated water level in February in Vejle Å. For the regularization term, we assumed that M of natural 473 rivers was characterized by spatial continuity and changes smoothly. However, most of the natural rivers have 474 been affected by human activities, for example, the construction of reservoirs and dams, vegetation management, 475 etc., which will often result in abrupt changes of hydraulic properties. Additionally, we introduced different 476 weights for regularization terms which are determined by trial-and-error method. A more comprehensive 477 method to select optimal weights for the regularization terms is the L-curve approach, which aims to 478 compromise data misfit and the regularization terms (Hansen and O'Leary, 1993). However, implementation of 479 such methods requires a large amount of computational resources. 480

Conclusion 481
This study presented a calibration approach for spatiotemporally distributed Manning-Strickler 482 coefficients M by combining UAS altimetry and in-situ stations. Several synthetic experiments under different 483 assumptions were evaluated and used to demonstrate the feasibility and efficiency of the scheme. We then 484 applied the approach in a real-world case study, an approximately 20 km river stretch of Vejle Å, to calibrate 485 and validate the M. 486 The synthetic study shows that M could be well-calibrated by the dense sampling of WSE in both spatial 487 and temporal scale. The synthetic calibration experiments showed that standard errors of ca. 3 cm on the WSE, 488 as achieved with our current UAS radar altimetry system, are sufficient to constrain spatio-temporal variations 489 of hydraulic roughness. Meanwhile, the steady-state solver is efficient and effective in simulating the WSE 490 along the river course, which is a requirement for the implementation of global search algorithms for highly The calibration algorithm was further applied on a real-world dataset to fit the optimal sets of M with 493 limited gauging data and UAS altimetry. The calibration results indicated that the M of Vejle Å changed 494 significantly in space and time. The simulated WSE has a RMSE of 8.48 cm on a high-flow day and 6.99 cm 495 on a low-flow day compared to highly resolved UAS altimetry. The simulated time series of WSE at four 496 hydrological stations have an average RMSE of 10 cm after we transferred the M calibrated by the steady-state 497 solver to the hydrodynamic model.

519
Scenario B1: Calibration of spatiotemporally varying Manning-Strickler coefficients M using in-situ 520 data only. This scenario is set for a river stretch monitored by discretely distributed monitoring stations. 521 Considering that the station data was continuously collected with a short time interval, lots of information will 522 be lost if we utilize only a few days for steady-state calibration. Here, the MIKE Hydro River model is used for 523 WSE simulating, and the local optimization is used for calibrating M. The initial value is essential for the 524 optimization, and the initial M was set to 20 m 1/3 /s for all the river stretches. Please be aware that it is extremely time-consuming if we set the calculation grid spacing as 20 m when using the MIKE Hydro River model for 526 Vejle case, thus propagations only calculated in each cross section. 527 Scenario B2: this scenario is the same as the former one, but with additional UAS altimetry collected 528 in February and June 2020 and partially covered Vejle Å, as shown in Figure 6b.