Estimating River Channel Bathymetry in Large Scale Flood Inundation Models

Flood inundation modeling across large data sparse areas has been increasing in recent years, driven by a desire to provide hazard information for a wider range of locations. The sophistication of these models has steadily advanced over the past decade due to improvements in remote sensing and modeling capability. There are now several global flood models (GFMs) that seek to simulate water surface dynamics across all rivers and floodplains regardless of data scarcity. However, flood models in data sparse areas lack river bathymetry because this cannot be observed remotely, meaning that a variety of methods for approximating river bathymetry have been developed from uniform flow or downstream hydraulic geometry theory. We argue that bathymetry estimation in these models should follow gradually varying flow theory to account for both uniform and nonuniform flows. We demonstrate that existing methods for bathymetry estimation in GFMs are only accurate for kinematic water surface profiles and are unable to simulate unbiased water surface profiles for reaches with diffusive or shallow water wave properties. The use of gradually varied flow theory to estimate bathymetry in a GFM reduced model error compared to a target water surface profile by 66% and eliminated bias due to backwater effects. For a large‐scale test case in Mozambique this reduced flood extents by 40% and floodplain storage by 79% at the 5 years return period. The wet bias associated with uniform flow derived channels could have significant implications for modeling the role floodplains play in attenuating river discharges, potentially overstating their role.


Introduction
In recent decades inundation modeling has become an integral component of flood management activities by providing hazard and risk mapping data to decision makers (Merz et al., 2010). Fundamental to this success has been the development of modeling frameworks where detailed river and floodplain bathymetry is combined with flow predictions or observations to construct numerical models that can simulate inundation depth for various scenarios (e.g., return periods) (de Moel et al., 2009). Results from such models have without question been useful to risk managers and are recognized by both national and international disaster risk reduction policies (Priest et al., 2016;Van Alphen et al., 2009). However, the expertize needed to implement these modeling methods, their data requirements and overall cost means that most flood risk data have been generated in developed countries, leading to substantial inequalities in risk information and an inadequate understanding of risk for many localities.
Perhaps unsurprisingly, there has been a move toward extending flood predictions to data and resource sparse areas, and in some cases using automated modeling approaches to enable regional or global coverage (Dottori et al., 2016;Paiva et al., 2011;Sampson et al., 2015;Siqueira et al., 2018;Ward et al., 2015;Winsemius et al., 2013). This has been supported by better numerical codes and improvements in the availability and accuracy of key datasets, such as global elevation models  and river width data (Allen & Pavelsky, 2018). Examples of regional and global scale inundation modeling now cover a range of applications including flood hazard estimation (Alfieri et al., 2014;Pappenberger et al., 2012;Sampson et al., 2015;Winsemius et al., 2013), flood event set and loss estimates (Quinn et al., 2019), flood risk and exposure modeling (Jongman et al., 2012;Ward et al., 2013), discharge estimation from remote sensing (Andreadis et al., 2007;Biancamaria et al., 2011;Durand et al., 2008;Neal et al., 2009), understanding wetland dynamics (Neal et al., 2012), estimating climate change impacts on flooding (Alfieri et al., 2017;Dottori et al., 2018;Hirabayashi et al., 2013) and modeling evaporative feedback to the atmosphere from wetlands (Dadson et al., 2010).
Despite wide ranging applications, the accuracy of global flood model predictions is poorly understood, but also likely to be low in many cases given that different models tend to disagree on where is at risk (Aerts et al., 2020;Bernhofen et al., 2018;Trigg et al., 2016), and this is particularly the case over complex floodplains such as deltas. Accuracy is known to be sensitive to DEM quality (Courty et al., 2019;McClean et al., 2020), vegetation and buildings removal from the DEM, extreme flow estimates, model resolution (Aerts et al., 2020;Fleischmann et al., 2019) and river channel cross-section accuracy (Fleischmann et al., 2019). Furthermore, GFM inundation outputs typically show less sensitivity to return period discharge than expected (Dottori et al., 2016) and tend to over-predict exposure for more frequent events (smaller magnitudes) relative to loss observations (Quinn et al., 2019), which is often assumed to result from a lack of flood defense information (Aerts et al., 2020;Ward et al., 2013). While not doubting the importance of flood defenses, and there are noteworthy efforts to improve these data (Scussolini et al., 2016), it is also imperative that the inundation model simulates an accurate water profile with respect to the river bank and floodplain heights. This is because most floodplains are inherently flat, and thus a small increase in simulated in-channel water height can generate a substantial increase in simulated flood extent. It follows that modest biases in the simulated water profile around bankfull discharge will adversely impact model accuracy during small high-frequency flood events. This has significant practical implications as high-frequency events must inherently make up the majority of events in any quantitative risk calculation (such as a loss-exceedance curve), and therefore they have the ability to significantly impact resultant risk estimates (e.g., Quinn et al., 2019). Furthermore, sensitivity analysis has shown that interactions between channel and floodplain due to uncertainties in channel bed elevation can influence wave propagation and flooding downstream (de Paiva et al., 2013).
This study will review the methods used for estimating river channel bathymetry when modeling floods in data sparse areas, focusing on their effectiveness at reproducing water surface profiles. Where the definition of the river channel typically includes properties such as channel location, bank full conveyance, bank height, width, depth and friction. Current methods are shown to be vulnerable to substantial errors when simulating river water surface profiles under nonuniform flow conditions, which we demonstrate results in an over-prediction of hazard. To address these issues and provide more robust simulation, improved approaches to river channel definition in the absence of cross-section data are proposed based on gradually varied flow theory. Gradually varied flow theory has been proposed for inverse problems in other hydraulic contexts, such as discharge estimation for the forthcoming Surface Water and Ocean Topography mission (e.g., Andreadis et al., 2020), and thus seem promising for improving flood predictions as well. Such approaches should provide more control over the behavior of channel-floodplain interaction in data sparse areas and lead to simulations that are more consistent with the wave theories that underpin the inundation modeling.

Approaches to Flood Inundation Modeling in Data Scarce Areas
In the most basic terms, all flood inundation models require at least four components, with the last item on this list being the focus of this study: 1. Boundary and initial condition that define the volumes of water flowing in the model domain. 2. A numerical modeling method to simulate river and floodplain flows. 3. Digital terrain data to define the floodplain surface. 4. A definition of the river network-specifically bank full conveyance and the subordinate variables of channel width, depth, section shape, and friction.
Each of these four components can be handled quite differently depending on the intended application of the model, the balance needed between compute speed and accuracy, the data available and the expertize of the model developers. For example, volume inputs can range from direct observation of past events at gauges (Pappenberger et al., 2006), to design hydrographs representing extreme events via a regionalization of gauging station data (Flood estimation handbook, 1999;Smith et al., 2015), to runoff inputs from hydrological and land surface models (Dottori et al., 2016;Paiva et al., 2011;Winsemius et al., 2013;Yamazaki et al., 2011). The numerical representation of river channel hydraulics in inundation models in global flood models also varies substantially in complexity, encompassing methods such as linear advection and diffusion wave methods (Lohmann et al., 1998), kinematic waves (Bell et al., 2007;Oki & Sud, 1998), diffusive waves (Sayama et al., 2012;Yamazaki et al., 2011), dynamic waves (diffusion + local inertia) (Neal et al., 2012;Yamazaki et al., 2013) and shallow water wave processes (Paiva et al., 2011;Sanders & Schubert, 2019). These generally increase in complexity and accuracy of process representation in the above order, with simpler methods generally applicable in fewer physical settings or used for large scale modeling.
To simulate floodplain inundation, the channel model must be linked to a model representing the floodplain conditioned on suitable digital elevation data (Courty et al., 2019;Ettritch et al., 2018;Hawker et al., 2018;Marks & Bates, 2000;Sanders, 2007;Sanders et al., 2005). Approaches to represent floodplain inundation at regional to global scale range in much the same way as those for the channel, from relatively simple DEM filling type methods (Nardi et al., 2019;Winsemius et al., 2013), where no dynamics are assumed, to extending one-dimensional (1D) model cross-sections onto the floodplain (UNISDR, 2015), to models that dynamically link the channel to large floodplain storage areas (Decharme et al., 2008;Paiva et al., 2011;Yamazaki et al., 2011) to models that simulate inundation dynamics in two-dimensions (Neal et al., 2012;Sanders & Schubert, 2019;Sayama et al., 2012). As with the channel models, complexity, cost and accuracy generally increase as you move down this list, with simpler models easier to apply over large areas. The spatial resolution of the simulations is usually governed by the resolution of the floodplain topography and a trade-off between acceptable computational cost and the spatial precision required by the application. Finally, the model will require a definition of the river channel network, which acts as a critical control on how water moves through a landscape. In the case of fluvial flooding, the river channel is usually the main conveyor of discharge and will interact with the floodplain in a complex manner as water moves both from and to the channel, depending primarily on topography and friction variability (Fewtrell et al., 2011;Knight & Shiono, 1996). Even in the case of other types of flooding, such as pluvial and coastal, the role of channel conveyance can be significant. How channels are represented in a data scare setting where bathymetry has not been observed will therefore influence inundation simulations significantly (Neal et al., 2012;Sampson et al., 2015;Yamazaki et al., 2011), particularly at low return periods where small changes in river conveyance can have a disproportionately large impact on the simulated flooding. For traditional reach scale hydrodynamic modeling, the quality of river bathymetry data, often in the form of channel cross-sections or surfaces from sonar data, is key to accurate simulation of the relationship between discharge and water level (Cook & Merwade, 2009). However, since such data are unavailable in data scarce contexts an approximation must be used that best represents the water surface elevation and discharge relationship given the available data (Grimaldi et al., 2018).

Methods for Defining River Channels in Data Sparse Flood Inundation Models
Several approaches have been proposed to simulating channel hydrodynamics in the absence of cross-section data. We review and categorize these starting from the simplest case of removing the channel component from the model entirely. Note that the methods cited also vary considerably in their treatment of hydrology, choice numerical scheme, and floodplain DEM, but these will be secondary considerations here.
For this discussion, we will take as a starting point that the profile of the river water surface can be defined as a gradually varied flow (GVF): where h is the depth, x is the distance downstream, Fr is the Froude number, 0 S is the bed slope and f S is the friction slope. If we also assume that the river channel is rectangular, and that friction is represented by Manning's equation then the friction slope is found via: where n is the Manning's roughness coefficient, w is the channel width and Q is the discharge. The Froude number of the channel Fr is then where g is the acceleration due to gravity. We argue that the inundation models need to accurately estimate the water surface profile p for bank full discharge Q bf along the river if the model is to simulate flooding during events or for specific design discharges. The simulation of the flow profile below bank full is not a priority in this case, justifying the commonly used simplification to a rectangular channel. Although, there will of course be trade-offs when seeking an accurate water surface profile at bank full that depend on the degree of data scarcity and if observations of the water surface and discharge are available below bank full discharge it might be better to estimate bathymetry from these rather than the bank full profile. However defined, obtaining an accurate flow profile from Equation 1 will depend on identifying the bed friction (n), bed elevation (z) from which bed slope ( 0 S ) is defined, channel width (w) and bank full discharge (Q bf ) along the river network. Developers and researchers have approached defining these in several ways as outlined below in order of increasing complexity.

The No Channel Method
The simplest approach to represent the river network is to estimate bank full discharge and then remove this from the event or design flood discharge of interest. This "excess discharge" is then used by a series of reach scale models of the floodplain without a river channel. This is a popular method due to the simplicity of not needing to estimate w, n, or z and can yield sensible results over large areas (Alfieri et al., 2014;Bradbrook et al., 2005;Dottori et al., 2016). However, floodplain flow pathways are usually complex and dominated by the interaction between the floodplain and channel (Lewin & Ashworth, 2014;Trigg et al., 2012). The nonlinear nature of the governing equations means that flood wave celerity will differ significantly if deep channel flows are omitted. Essentially fast-moving water in the channel interacting with slow moving water on the floodplain is needed to correctly simulate wave propagation during a flood. This approach is also particularly sensitive to the presence of objects or errors in the topography that impede the propagation of the flood wave downstream because the channel network would mitigate for these by allowing water to return to the channel and move on downstream at greater velocity. Over prediction as a result of mass blockage effects is therefore a concern (e.g., Neal et al., 2012) such that this approach works best for simulating short reaches using very accurate terrain data for example, the UK extreme flood zone maps of Bradbrook et al. (2005). Over large and complex floodplains and deltas, we expect the method to become inaccurate (Neal et al., 2012;Sampson et al., 2015). For example, Fleischmann et al. (2016) shows how channel floodplain interactions on the Amazon river are necessary to recreate the skew and arrival time observed in downstream gauge data.

Empirical/Hydraulic Geometry Methods
An alternative to removing the channel is to estimate its dimensions empirically given observations from surveyed rivers using downstream hydraulic geometry theory, as demonstrated in the GFM of Yamazaki et al. (2013) and many regional models (Fleischmann et al., 2019;Grimaldi et al., 2019;Mateo et al., 2017;Neal et al., 2012;Paiva et al., 2011;Sayama et al., 2017;Schumann et al., 2013). Downstream hydraulic geometry theory aims to estimate how the width and depth of the channel are related to bank full discharge by a series of power laws (Leopold & Maddock, 1953). These parameters have been estimated empirically from field observations over many sites to allow prediction at locations without observations Hey & Thorne, 1986). Once depth has been established, the bed elevation is typically calculated by subtracting depth from the river bank height all along the river network as defined in the DEM (with the aid of some processing along the channel to reduce DEM noise Yamazaki et al., 2013Yamazaki et al., , 2019. This method benefits from being simple to implement but there is little chance that the desired surface profile will be simulated at bank full discharge because, (a) the hydraulic geometry parameters are uncertain and difficult to regionalize, with substantial variability expected between rivers; (b) the friction value in the hydrodynamic model will need to be estimated because there is no direct link between the hydraulic geometry parameters and friction parameters; and (c) changes in profile slope are not accounted for, meaning the hydrodynamic model will simulate a different water surface elevations to those expected.
Given that width is readily observable from remote sensing platforms (Isikdogan et al., 2017;Lin et al., 2020;Yang et al., 2020) versions of this approach where the widths and other observable factors are used to help predict the bed elevation or bank full discharge have also been proposed (Gleason & Smith, 2014;Grimaldi et al., 2018). This approach has the advantage of not needing an estimate of bank full discharge, however when width is used to predict the depth a reach will shallow when the channel narrows and deepen when it widens if widths are not appropriately reach averaged, thus changing the conveyance in an unrealistic manner. This method was used by Neal et al. (2012) over a delta where the proportions of flow bifurcating down tributaries were unknown and strong evaporative feedback meant that mass was not conserved along reaches. However, the power law parameters and model friction were so uncertain that they needed to be estimated from water surface observations via a computationally expensive calibration process.

Uniform Flow Theory-Manning's Equation Method
A simple way to calculate the channel bathymetry is to assume that a uniform channel exists over long distances when calculating the depth such that uniform flow formula can be used (e.g., Miguez-Macho & Fan, 2012;Sampson et al., 2015). Under this assumption the bed slope 0 S and friction slope f S are assumed equal, and if the channel is further assumed to be sufficiently wide that hydraulic radius (cross sectional area divided by wetted perimeter) is equal to depth (cross sectional area divided by width) then the channel depth h can be calculated analytically using Manning's equation from the local water surface slope S, friction n, width w, and discharge Q.
Following the GFM described by Sampson et al. (2015) the channel bed is found by subtracting depths from a smoothed bank high profile. This approach overcomes some of the limitations of the hydraulic geometry method by accounting for the friction parameterization of the model (i.e., the friction might not be known but at least the same value can be used for the flood model and bathymetry estimation) and allowing depth and width to vary inversely for the same bank full discharge. However, if observations of the water surface heights are used to define the slope S the discharge must be known at the same time as the observations were taken, while in the case where bank heights are used to derive slope an estimate of the bank full discharge is needed as with the hydraulic geometry methods. Furthermore, we know that many controls exist on flows that cause the river to depart from uniform conditions, for example constrictions in channel width, changes in discharge (e.g., a tributary joining), changes in bed slope and the presence of water bodies such as lakes. Backwater effects from these controls will be significant in most lowland channels (Trigg et al., 2009) meaning that the uniform flow assumption in Sampson et al. (2015) GFM could underpredict the size of the channel needed to simulate the correct water surface at bank full discharge during hydrodynamic simulation.

Nonuniform Flow
A more accurate method than those above would be to use a bed profile that simulates the desired water surface profile given the gradually varied flow equations themselves, such that nonuniform flow is accounted for. Unlike the uniform flow case there is no analytical solution, however solving these equations is well established when river bathymetry is defined (Chaudhry, 2008). The calculation process involves starting from a control section where the water level is also known and then integrating upstream in the case of subcritical flows or downstream in the case of supercritical flows (Henderson, 1966). In the GFM case, the profile has been observed rather than the bathymetry and the control section can be a lake, the ocean or river where uniform flow is assumed. To estimate a channel for a nonuniform flow profile we must therefore find the channel bathymetry that best approximates the observed profile given Equation 1. Several methods have been proposed to achieve this aim given data from satellite altimetry, shoreline observations and the proposed Surface Water and Ocean Topography mission (Andreadis et al., 2020; Brêda et al., 2019; Durand et al., 2008Durand et al., , 2014Garambois & Monnier, 2015;Yoon et al., 2012) which are covered in more detail in Section 3.1.1.
The next sections present two bed estimation methods targeted at two test cases. The first test evaluates the principal of using nonuniform rather than uniform flow theory to estimate the channel bed, using a very simple, but well studied, real world test. It considers a reach scale situation where the water surface profile has been observed at a specific time and the discharge is known. The purpose of this test is to compare with observed bathymetry and assess the potential for improved accuracy at reach scale. The second test demonstrates practical application to a GFM. For this test we implement the bed estimation globally by simplifying the approach in test 1. The GVF based method is then benchmarked against the Manning's equation method of Sampson et al. (2015) for a test case in Mozambique and Malawi.

Test Case 1: Reach Scale Bed Estimation From Water Surface Elevation
The data for this test case were obtained for a widely studied reach of the River Severn (UK) that flows from Worcester to Tewkesbury (Bates et al., 2006;García-Pintado et al., 2015;Neal et al., 2015;Schumann et al., 2009). An observed water surface profile was sampled at 10 m intervals from a 0.5 m resolution airborne LiDAR survey conducted by the Environment Agency on the December 12, 2014 (similar to the approach of Smart et al. (2009)). LiDAR water surface returns for this reach are expected to include vertical error of 5-15 cm, equivalent to another flat surfaces. River discharge at the time of the LiDAR acquisition was measured at the Saxons Lode gauging station, around the middle of the study reach (See Figure 1 for map). Flows were within bank and assumed to be without error for the purpose of this test given the high-quality ultrasonic gauging station installed at this site and its use for operation flood forecasting (Q error <10% given the flow conditions). Ideally the water surface profile would be as close to bank full as possible to minimize uncertainties from the channel geometry between the observed profile and bank full profile when the channel geometry is used for flood modeling, however, in practice the timing of the LiDAR survey will be based on other priorities. The gauged discharge of 225.6 m 3 s −1 was exceeded only 5% of the time at this station and was assumed to be constant along the 25 km reach at the time of the LiDAR survey. The average water surface slope was 0.0001 mm −1 . Observations of the bed elevation and water surface top width were obtained from cross section data provided by the Environment Agency. Manning's roughness coefficient was assumed to be 0.035 unless otherwise stated, which is physically reasonable for a reach like this with a gravel bed and cohesive banks. As flows in this reach are always subcritical, the GVF solver requires a water surface elevation and depth at the downstream boundary. These were estimated from the LiDAR observations of water surface elevation, slope and channel top width using Manning's equation.

Estimating Bathymetry
The bed estimation method presented below derives from recent work, in anticipation of the NASA Surface Water and Ocean Topography (SWOT) mission, which has focused on joint estimation of Q, n, and z from a time series of water surface height and slope observations (Andreadis et al., 2020;Durand et al., 2014Durand et al., , 2016. Specifically, the method presented here is a simplification of that proposed by Garambois and Monnier (2015) for SWOT discharge and bathymetry estimation, which aimed aim to estimate Q, n, and z from multiple observations of water surface elevation and slope through time. In our implementation Q and n will be known to the algorithm and there is no time element to our approach. Thus, our approach assumes discharge varies in space, but that the river can be approximated as steady state. This steady state assumption could be relaxed, at considerable computational expense, with the use of a 1D hydrodynamic model and time varying discharge in place of the gradually varied flow solver. Brêda et al. (2019) evaluate several data assimilation methods that would be suitable for bed estimation in such circumstances, thus these are not discussed here. Meanwhile there are a range of alternative methods for finding the inverse of the gradually varied flow equations or similar, which are reviewed in detail by Sellier (2016) but not tested here. We also assume that width can be measured accurately from the cross-section data and that the rectangular channel assumption is a reasonable approximation of the channel shape in this case.
The steps taken to estimate the channel bed are outlined in Figure 2 and in detail below. A first-order estimate of the bed elevations are needed as initial conditions, which would most obviously come from the Manning's equation method (uniform flow assumption) described above. From these bed elevations the gradually varied flow profile given the Manning's method bed is found using the Runge-Kutta method to solve Equation 1 (function ode45 in Matlab).
To refine the estimate of river bed elevations z from the first-order approximation, we seek the bed elevation that minimizes the least squares difference between the desired water surface profile p and that simulated by the gradually varied flow solver. The desired profile would ideally be from water surface observations, as in the Severn test here, but might also refer to bank or defense heights when implemented in a data scarce region. The response to changing bed elevation will be complicated due to backwater effects, therefore, nonlinear least squares optimization was undertaken using the trust region reflective method (Moré & Sorensen, 1983) to search within pre-defined bounds (z ub and z lb ) for the bed elevation, where subscripts ub and lb signify the upper and lower bounds respectively. The optimization function used here (lsqnonlin in Matlab) requires a vector of residual values (differences between observed and simulated profile) as input from which to compute the sum of square errors, therefore for the case where the water surface profile p is being estimated given discharge Q the following vector values are to be optimized NEAL ET AL.
10.1029/2020WR028301 8 of 22 Figure 2. Conceptual diagram of a bed estimation method for test 1, using a gradually varied flow solver and non-linear least squares estimation approach.
where θ contains all other parameters of the gradually varied flow solver for example, roughness and channel widths and   , , M z Q θ is the water surface elevation simulated by the gradually varied flow solver. This condition was enough to find a channel bed that simulates a water surface profile in this test. However, the solution is non-unique and in practical applications this method alone created undesirable outcomes in two key ways. First, the bed elevations varied from one estimation point to the next far more than expected, leading to a bed that was clearly physically unrealistic (discussed later in Section 3.1.2 and Figure 3). Essentially the magnitude of bed slopes is greater than expected and often includes substantial adverse slopes that are implausibly steep and inconsistent with natural bedforms. The water surface slope could also change substantially from observation to observation in an unphysical manner, suggesting an overfitting to the noise in individual water surface observations. Second, the optimizer failed to find the minimum root NEAL ET AL.
10.1029/2020WR028301 9 of 22 Optimized bed elevations and simulated water surfaces relative to LiDAR observations for (a) GFV solver method with no regularization terms, (b) GFV solver with shallow bed preference, (c) GVF solver with low gradient bed preference, and (d) GVF solver with low gradient water surface constraint. Also plotted are the initial bed estimated from Manning's equation and the simulated water surface from this. The channel thalweg (lowest point of channel bed) is plotted as crosses where available and the observed channel average bed is plotted as points given the hydraulic depth (cross sectional area divided by top width) subtracted from the LiDAR height. GVF, gradually varied flow. mean squared error (RMSE) (discussed later in Section 3.1.2 and Figure 4), indicating that it was unable to find a global minimum.
To provide a physically plausible bed and easier to optimize solution we evaluate the effectiveness of including three additional costs to the objective function to regularize (simplify) the solution, which aim to reduce the variability in bed elevation from one location to the next by penalizing for greater bed gradient (z), water surface gradient (     z h ) or channel depth (h). In other words, we are adding regularization terms to the optimization to prefer simpler solutions that are shallow or give gradually varying channel beds or water surfaces by passing one of the following vectors to the optimizer: How much weight is given to each of the additional costs is determined by parameter ρ, where the greater the parameter value the greater the weight and greater preference for a simple solution to the bed. The weight given to the gradient based costs needs to vary depending on the slope of a reach S to have a consist-NEAL ET AL.
10.1029/2020WR028301 10 of 22 ent effect against any given error magnitude. For example, we expect lower gradient rivers to require greater weights. Therefore, we normalize for this effect by multiplying ρ by the inverse of river slope S in these cases. S can be estimated from the longitudinal water surface, bank or floodplain elevations. In the next section, we will test this gradually varied flow solver approach and examine the effect of ρ for cost functions in Equations 6a-6c on the optimization routine.

Results From Test 1
A series of tests were implemented to assess different configurations of the bed solver objective function and establish the performance of the method verse assuming uniform flow. The simplest configuration seeks the minimum between the simulated water surface and the LiDAR water surface, while the others give differing weights to the preference for a low gradient and shallow bed and a low gradient water surface as described by Equations 6a-6c. Upper and lower bounds for the bed (z ub & z lb ) were set such that the channel could be no more than 30 m deep from the observed water surface and no higher than 5 m above the observed water surface. This range goes well beyond physically plausible values for depth for this reach which is typically in the range 5-8 m. Manning's equation (Equation 4) was used to estimate the initial bed profile and to provide a benchmark. Details of the Manning's equation implementation are provided in the supplement, however when the Manning's bed profile was used as input to the gradually varied flow solver, as defined in Equations 1-3, a RMSE of 0.34 m and mean error (ME) of −0.24 m to the LiDAR water surface was achieved.
Estimating the bed elevation from the LiDAR water surface, using Equation 5, produced a water surface with a RMSE of 0.16 m and ME of 0.073 to the LiDAR observations, 53% more accurate and significantly less biased than the Manning's based approach. The bed profile is plotted in Figure 3a and shows a high degree of bed variability relative to the observed bed data, with many locations where the channel is unrealistically deep relative to the observed bed. The observed bed data is presented in two way in Figure 3. First, crosses denote the lowest point of each irregular section, which we expect to be below the GVF solver estimates because a rectangular channel has been assumed. Second, the average observed bed was calculated given the observed hydraulic depth and plotted as dots to facilitate direct comparison with the GVF solver bed estimates. The hydraulic depth is the irregular section area below the LiDAR divided by top width used by the GVF solver.
Figures 3b-3d plot bed estimates and corresponding simulated water surfaces for the shallow channel (3b) bed gradient (3c) and water surface gradient (3d) based regularization approaches, given a value of ρ that gave the lowest RMSE. All the regularization approached resulted in a simplification of the bed elevations relative to the no regularization case in Figure 3a, albeit the impact being minor in the case of the low gradient water surface constraint. Therefore, the regularization was successful in reducing the complexity of the riverbed profile and made it easier to find a solution to the observed water surface.
A comparison of the average observed bed at each cross-section with the estimated values confirms that the simplification of the bed profile was also a more accurate approximation of the observed bed. Table 1 shows the root mean squared error (RMSE) and mean error (ME) for each of the bed profiles plotted in Figure 3, with the preference for a low gradient bed performing best (RMSE 0.77 m). Much of the RMSE derives from bias, which indicates a general overprediction of the bed in the case of the shallow bed and low gradient bed preference (0.37 and 0.38 m respectively) and a more substantial under prediction by the other solvers (−1.73 and −1.07 m). For the shallow bed and low gradient bed estimates increasing the friction parameter could NEAL ET AL. All statistics relate to the model data plotted in Figure 3. GVF, gradually varied flow.

Table 1 Accuracy of GFV Solver Bed Estimates Given Different Regularization Schemes
compensate for the mean error and further reduce the RMSE. However, without the regularization it was not possible to avoid very low bed estimates and a general underprediction of the bed elevation occurs at all plausible friction values (not plotted).
The impacts of considering bed slope, water surface slope or bed depth in the cost function (Equations 6a-6c) are further plotted in Figure 4, where (4a) plots RMSE against parameter ρ for each regularization term, (4b) plots the mean errors and (4c) plots the number of iterations the nonlinear least squares optimizer needed to find a solution, effectively the relative computational cost. The most accurate of the regularization terms preferred a lower gradient bed, which could reduce the RMSE to 0.095 m and the ME to −0.007 m for a ρ value of 0.046. This is a 72% reduction in RMSE relative to the Manning's model and is within the vertical error of the LiDAR survey. The fact that when a regularization term is included in the optimization the RMSE relative to the LiDAR water surface is reduced (0.16-0.095 m) indicates that without regularization, the optimizer fails to find a global optimum given the uniform flow based initial bed estimate. Regularization thus appears to help the numerical solver, in addition to providing a more physically realistic solution.
Imposing a cost for deeper beds reduced the RMSE to 0.099 m with a mean error of 0.014 m for a slightly lower ρ value of 0.033, which at first glance suggests this regularization approach could be almost as successful as the constraint on bed gradient. However, the range of ρ that produced good results (defined as results that are better than the approach without the regularization) is much narrower and mean errors (over-prediction of the water surface) increase with the magnitude of ρ. Regularizing on the water surface gradient could also reduce RMSE with higher values of ρ, however not by as much as the other two methods. This approach also required over twice the computational cost at the respective optimal values of ρ, suggesting it was better to regularize directly on the variable being estimated (channel bed elevations). Overall, the most accurate regularization approach in this case was to include a preference for a low gradient bed, which reduced water surface RMSE at no additional computational cost compared to not including the regularization terms. Values for ρ > 0 and <0.17 would all improve the simulation accuracy and bias, however the greatest reduction in RMSE occurred for a value of 0.046, indicating that scaling the bed gradient to around half the observation errors in the water surface profile produced best results in this case. There might be benefits in combining the constraints, this is possible, and we did try it. However, RMSE was not improved with more complex regularization schemes, which was expected given the observation errors, while the inclusion of the shallow bed preference always increased the mean error due to the shallow bias. Further reach scale tests for a range of flow profiles would be needed to test the robustness of this conclusion.
The introduction of regularization trades off minimizing RMSE for the other desirable properties (a smooth bed etc.) meaning the optimal RMSE is not found. It is likely that RMSE can be reduced further by reinitializing the solver with ρ = 0 using the regularized bed solution as a starting point, and this might find a global minimum RMSE. However, in this case the RMSE is already around the error magnitude we expect in the LiDAR observations and in our opinion, there is good argument for accepting a sub optimal RMSE to maintain the other desirable properties of the solution and mitigate for overfitting to the profile data.
Finally, additional tests were conducted given channel friction values from 0.02 to 0.07 s/(m 1/3 ), a range of initial starting depths from 1 to 20 m, and a range of model resolutions from 100 m to 2,000 m. In all cases the bed gradient regularization term ρ was set at 0.046 (the optimum from the previous experiments). The results of these experiments are tabulated in Table 2. Friction has a negligible impact on RMSE and ME, although lower friction values required more function calls to the GVF solver in order to perform the optimization, potentially because the channel is shallower and further from the first-order estimate NEAL ET AL. Function calls to the GFV solver from the least squares nonlinear optimizer are also shown and are proportional to computation cost (note that the higher resolution models are also more expensive due to the number of estimation points). GVF, gradually varied flow; ME, mean error; RMSE, root mean squared error.

Table 2 Impact of Changing Channel Friction, Initial Depth of Channel and Resolution on GVF Solver Water Surface Accuracy (In Terms of RMSE and ME to LiDAR Water Surface Observations)
based on Manning's equation. The initial bed elevation had no impact on RMSE or ME, demonstrating that the solver did not get stuck in local optima for this test. A better initial bed required substantially fewer function calls (around half) of a poor initial bed (e.g., 1 m or 20 m), but if Manning's equation is used for the initial bed and the friction is the same between the two method it is unlikely there will be many cases where the initial bed is as poor as some of those tested here. The resolution of the observations and bed estimation points influenced the accuracy of the results with respect to both RMSE and ME, with higher resolutions being more accurate as expected. For this reach, RMSE increased by a little over 0.02 m from 2 to 1 km, but by only 0.1 cm from 1 km to 100 m. The number of function calls increased at finer resolution due to the need to estimate more bed elevations, between 2 and 0.5 km the function calls per estimation point were essentially the same, but doubled for the 0.1 km resolution indicating that it was more difficult to fit finer resolution data.

Test Case 2: Implementation in a GFM
Given the improved profile accuracy seen at the reach scale, a bed estimation method based on gradually varied flow was implemented within a GFM Wing et al., 2017). This model previously used the Manning's method for bed estimation and is described in detail by Sampson et al. (2015). A full description of the GFM is beyond the scope of this paper; however the key components include: (a) a regional flood frequency analysis  to provide return period discharge for all points on the global river network; (b) river network and terrain data sets based on the MERIT DEM  and MERIT HYDRO (Yamazaki et al., 2019) from which river locations and floodplain elevations were extracted; and (c) a regionalized river width estimation approach as described by Sampson et al. (2015). Numerical simulations were performed using the LISFLOOD-FP hydrodynamic model given the 1D channel from the bed solver and a two-dimensional (2D) floodplain model based on MERIT DEM. LISFLOOD-FP is a hydrodynamic model that solves a simplification of the shallow water equation without convective acceleration terms (Bates et al., 2010;de Almeida et al., 2012). For large scale applications, it uses a regular gird with either geographic (WGS84) coordinate systems and a 1D sub-grid scheme for river channels (Neal et al., 2012;Sampson et al., 2015).
The GFM was implemented at 3 arc second (∼90 m resolution) with bed estimates linearly interpolated from a 30 arc second (∼900 m resolution) river network that includes all rivers with an upstream catchment area above 50 km 2 . In the previous test case the profile p, which the channel bed was optimized to, was defined from observations of water surface elevation at a particular time associated with a gauged flow. In this "data sparse" test case the profile must be defined from bank heights and associated with an estimate of bank full discharge based on return period because instantaneous water surface elevations and discharges are unavailable. For small channels, the initial bank elevation profile p was taken from DEM cells directly above the channel and conditioned using a local smoothing function and monotonicity constraint. For larger rivers, a sample of elevation values for each node was taken from DEM cells adjacent to the water mask (SRTM water body data) to ensure that bank elevations (rather than water elevations) were being sampled. A clustering algorithm was used to segment the sampled elevations, with the median value of the lowest cluster assumed to be the bank elevation. As with small rivers, the channel elevation profile was then conditioned using a local smoothing function and monotonicity constraint . For some rivers, returns from sand bars and other feature below bank full will be included in the bank high estimates, where these are infrequent they will be ignored by the median filter, however it is likely that in some places the bank heights are reduced by these data. All rivers were assumed to convey the 1 in 2 years flow at bank full discharge, which is typical of most GFM implementations. Flood inundation extents are expected to be sensitive to this assumption, thus, although we aim to improve the simulation of the bank full water surface profile by the GFM for a given return period, we do not reduce the uncertainty in the bank full discharge, which can be substantial (Williams, 1978). Bank full return period might be refined using data on historical flood occurrence  regionalization of section observations (Castro & Jackson, 2001) or using standard of protection for engineered reaches (Bates et al., 2021;Winsemius et al., 2013).
For global application there are ∼0.5 billion locations where the channel bed must be estimated. Thus, although it is in theory possible to implement the GVF method as applied to the Severn test case globally by splitting the worlds river network into many small reaches and using downstream reaches as boundary conditions for those upstream, computational constraints required a simplification to the estimation process whereby the trust region reflective nonlinear least squares estimation was replaced by a simpler bed nudging approach to optimize the bed elevations from the initial bed. This nudging involved the following steps: i) Solve for the water surface profile given the initial bed estimate from Manning's equation ii) Compute the differences between bank profile and simulated water surface elevations from step i iii) Apply the differences between the bank profile and simulated water surface from ii to the bed elevations iv) Recompute the water surface profile, and v) Repeat ii, iii, iv once more to get a final set of bed elevation estimates and water surface profile errors.
We expect this method to be less accurate than the global optimization approach taken on the Severn. For example, our nudging approach is similar to the "direct insertion" method proposed by (Brêda et al., 2019), which they found to be less accurate that global optimization and Kalman filter based techniques technique when assimilating altimeter data.
The GFM is implemented as 10° × 10° overlapping tiles. Downstream boundary conditions for each river were estimated using Manning's equation at the edge of each model tile or using the GVF method results for the river mainstem where tributaries join. Coastal water heights were set to mean sea level. For sections of supercritical flow, the channel depth was set to the critical depth because the same Froude limit is applied in the LISFLOOD-FP inundation model when implemented at large scale following Adams et al. (2017). A gradually varying flow solver capable of simulating supercritical flow transitions could be used with shallow water hydrodynamic models. A convenient implication of this assumption is that the GVF solver only ever needs downstream boundary conditions as the supercritical profile is never simulated.
One GFM tile has been chosen here for further analysis to keep the data volume presentable in the following plots. It covers 10°-20° south and 20°-30° east, which includes North and Central Mozambique, Malawi, southern Tanzania and the eastern edges of Zambia and Zimbabwe. Inundation data are presented for a smaller 3° × 3° region including the town of Beira and the mouth of the Zambezi river for visualization purposes ( Figure 5). Beira experienced substantial flooding in 2019 due to cyclone Idai, which was the deadliest cyclone on record to have hit Africa. The GFM data were used to produce disaster bulletins that include estimates of population exposure as outlined by Emerton et al. (2020). The return period of the flooding is unknown, however satellite imagery of the flooding from the Copernicus Emergency Management Service (sentinel 1 data) best fits the GFM inundation extents at the 250-years return period. Thus, the GFM simulations are consistent with the meteorology of this rare event. Since the site includes the delta of the Zambezi, where the channel can be over 1 km wide, headwaters and extensive floodplains/wetlands we believe it is indicative of most locations where GFM data may be used.
Results from the Manning's and GVF method are presented in Figure 6, which plots the difference between the water surface profiles at bank full discharge simulated by the two methods and bank height for all 144,523 channel locations within the domain. This accuracy measure is plotted against a) bank height, b) bankfull discharge, c) bank slope, and d) Froude number (defined by Equation 3) in order to understand how the wave properties and physical setting affect water surface profile accuracy.
Overall, the root mean squared error between simulated water surface elevations and bank heights was 0.872 m for the Manning's method and 0.291 m for the GVF method. Mean error was 0.167 m for the Manning's method and −0.030 m for the GVF method. From these numbers and visually in Figure 6 it was clear that the Manning's method tended to over predict the water surface relative to the target bank heights, whereas the GVF solver was unbiased while also more accurate. Elevation and bank full discharge had no systematic impact on the magnitude of errors for either method. However, Manning's method errors were generally larger for slopes less than 10 −3 m m −1 , with almost all of the overprediction at low Froude (<0.2). The GVF method was almost always more accurate than the Manning's method, however there were 16 points with Froude number >0.99 where errors exceed −2 m. These points of poor performance correspond to steep features such as dam outlets and waterfalls, however the impact was limited locally to just a few model cells and did not propagate widely through the model domain. These errors are likely a weakness of out GVF solver setup in that we do not consider supercritical flow profiles or discontinuities.
Given the poor performance of the Manning's method at low Froude we investigated the relationship between water surface error and wave type following the approach of Trigg et al. (2009) following Vieira (1983). This was done by calculating the kinematic wave number k in addition to the Froude number at every river location: NEAL ET AL.
10.1029/2020WR028301 15 of 22 where L is the channel length (1 km in this case). When plotted against Froude number, suitable approximations for the wave at each point on the river network can be identified. According to the analysis of Vieira (1983), a kinematic wave is considered to be a reasonable approximation for locations with approximately k > 10 and Fr > 0.5, for k > 3 and Fr < 0.5 the wave can be characterized as diffusive, while shallow water characteristics are important for lower values of k. Figure 7 plots each river location into this wave characteristics space for both bed estimation methods. The dot color indicates the magnitude of the profile error (truncated at 1 m to aid visualization). Manning's method results (Figure 7a) were only accurate for reaches that can be approximated by a kinematic wave, which was expected because uniform flow is assumed. The channel depth was underestimated where diffusive or shallow water wave processes become important. The poor performance of kinematic wave models over large lowland rivers and the important role backwater effects can have on flooding is well established (Ikeuchi et al., 2015;Trigg et al., 2009;Yamazaki et al., 2011) and these results are in line with the expectation that the Manning's bed estimation method would over-predict the water surface and flood extent.
The GVF method was more accurate when diffusive and shallow water wave characteristics were important, with errors >1 m mainly present at the very lowest Froude and kinematic wave numbers. Further investigation revealed these locations to be lakes (predominantly lake Malawi). It is important to note that the GFM does not include a bespoke routine for simulating lake levels. Rather, lakes are represented as very flat rivers with widths as defined as if they were rivers . Given this combination of unrealistic lake widths and a flat bank it was not possible to accurately simulate the lake surface profile, and in classic applications of GVF theory lakes would usually form boundary conditions rather than be estimated. Furthermore, tributaries joining the lake connect to the flat mainstem moving down the center of the lake account for the encircled line of negative bias points (∼−1 m) at a range of different slopes in Figure 6c.
NEAL ET AL.
10.1029/2020WR028301 16 of 22 The range of channel slopes for the lake are due to the transition from the steep slopes either side of Lake Malawi to the flat mainstem in the middle of the lake. Results for lakes are usually discarded and replaced by a water mask when postprocessing the GFM data, meaning the impact on risk estimates from the GFM is likely to be small. However, an improved consideration of lakes will be needed for the GFM to simulate lake levels. There are also some cases of positive errors in the shallow water flow zone, these indicate the potential for reduced accuracy under such conditions. However, the same errors are present in the Manning's method used for initial conditions, which are not corrected by the simple optimization routine used here. Implementing the nonlinear least squares estimation method from test 1 or supercritical flow profiles may be necessary in some location to gain accurate results.
Although errors in the water surface profile can be several meters greater for the Manning's method than the GVF method we do not expect such large errors in water surface elevation to propagate into the inundation extents because floodplain storage and conveyance will dampen the wave amplitude. Figure 5 plots the flood inundation extents for the 1 in 5 years return period flood to assess the impact of the surface profile errors on inundation simulation. Flood extent was 40.3% greater for the Manning's method and total floodplain storage increased by 79.4%. Meanwhile exposure to this return period was 150,000 people for the Manning's method and 90,000 people for the GVF method according to High Resolution Settlement Layer (HRSL) population data (https://www.ciesin.columbia.edu/data/hrsl/). Therefore, inaccurate specification of the channel bed by not accounting for nonuniform flow will bias flood inundation and exposure calculations. Studies that use such methods for management activities such as assessing the value of floodplains for wave attenuation and national/international scale overviews of flood defense requirements (especially for low return period events) should take note of these substantial biases. For greater flows, the subtle differences in bed profile are substantially less important and changes in inundation extent are locally complex. For the 100-years flood the Manning's method increases the volume on the floodplain by 5% but decreases the flood extent by a similar percentage. Population exposure to this return period was 431,000 people for the Manning's method and 475,000 people for the GVF method because more of the water tends to go out of bank further downstream in the GFV channel. Validation of these extreme inundation extents is beyond the scope of this paper, however an extensive analysis of the GFM hazard estimates using this solver is presented by Bates et al. (2021).
NEAL ET AL.

Conclusions
This paper has developed and demonstrated methods for channel bed estimation based on a simple gradually varied flow solver, which are suitable for application in reach scale and global scale flood models. The principal of considering nonuniform flow rather than uniform flow when estimating the bed elevations was first evaluated on a well-studied site in the UK where the accuracy of the approach could be assessed, and different numerical methods evaluated. We found that: 1. The GFV method outperformed a Manning's equation method and could reconstruct a LiDAR observed water surface profile to within the expected observation errors (<10 cm RMSE). 2. Regularization was necessary to provide realistic bed profiles. For our test, an effective way to do this was to add a cost to the objective function for stepwise changes in bed elevation scaled in proportion to the reach slope. 3. The method was robust when given friction parameterizations within typical ranges and poor initial bed estimates. 4. The bed estimation process also performed well across resolutions from 100 m to 2 km, although cost increased rapidly for negligible accuracy gain toward 100 m resolution and accuracy was reduced between 1 and 2 km resolution.
A simplified GVF bed estimation method was implemented in a global flood model .
Results from a test case in east Africa demonstrated that model errors to a defined water surface profile were reduced from 0.872 m for the Manning's method to 0.291 m for the GVF method. Bias toward over prediction by the Manning's method was also eliminated with mean error falling from 0.167 m for the Manning's model to −0.030 m for the GVF method. Improvements over the Manning's method were greatest for profiles with diffusive or shallow water wave properties, highlighting the importance of backwater effects on the flow profile.
The improved bed estimates had a substantial impact on floodplain inundation and storage dynamics for small floods when use in a hydrodynamic model. For the 1 in 5 years return period, inundation extent was reduced by 40.3%, total floodplain storage decreased by 79.4% and exposure fell from 150,000 people to 90,000 people over the study domain. Therefore, inaccurate specification of the channel bed by not accounting for nonuniform flows biased flood inundation, storage and exposure calculations. Although bed estimates were less important for larger magnitude flows.
Studies that use such methods for management activities such as assessing the value of floodplains for wave attenuation and national/international scale overviews of flood defense requirements (especially for low return period events) should take note of these biases as they might be significant (not tested here). Similarly, any flood defenses added to the model will likely underperform during hydrodynamic simulation for the same reasons. For financial services applications, where catastrophe risk modeling based on event sets is often needed, the issue is particularly acute. This is because event sets must simulate flooding from both large and small return period flows in order to estimate key risk metrics such as loss exceedance probabilities. These are unlikely to be accurate with channel geometries based on Manning's equation or other simpler methods because of too much flooding for lower return period flows.