Efficient incorporation of channel cross-section geometry uncertainty into regional and global scale flood inundation models.

This paper investigates the challenge of representing structural differences in river channel cross-section geometry for regional to global scale river hydraulic models and the effect this can have on simulations of wave dynamics. Classically, channel geometry is deﬁned using data, yet at larger scales the necessary information and model structures do not exist to take this approach. We therefore propose a fundamen-tally different approach where the structural uncertainty in channel geometry is represented using a sim- ple parameterisation, which could then be estimated through calibration or data assimilation. This paper ﬁrst outlines the development of a computationally efﬁcient numerical scheme to represent generalised channel shapes using a single parameter, which is then validated using a simple straight channel test case and shown to predict wetted perimeter to within 2% for the channels tested. An application to the River Severn, UK is also presented, along with an analysis of model sensitivity to channel shape, depth and friction. The channel shape parameter was shown to improve model simulations of river level, particularly for more physically plausible channel roughness and depth parameter ranges. Calibrating channel Manning’s coefﬁcient in a rectangular channel provided similar water level simulation accuracy in terms of Nash–Sutcliffe efﬁciency to a model where friction and shape or depth were calibrated. However, the calibrated Manning coefﬁcient in the rectangular channel model was (cid:2) 2/3 greater than the likely phys- ically realistic value for this reach and this erroneously slowed wave propagation times through the reach by several hours. Therefore, for large scale models applied in data sparse areas, calibrating channel depth and/or shape may be preferable to assuming a rectangular geometry and calibrating friction alone. (cid:2) 2015 The Authors. Published by Elsevier B.V. ThisisanopenaccessarticleundertheCCBYlicense(http:// creativecommons.org/licenses/by/4.0/).


Introduction
Recently there has been substantial interest in simulating river hydraulics at regional to global scales, most notably for the purpose of flood hazard and risk assessment. There is currently no clear definition of what constitutes a large scale hydraulic model, however for the purpose of this paper it will be assumed that the model has a structure that can be applied to simulate water levels and flows over an entire continent. Implicit to this definition is that some form of specialisation of the model structure will have occurred to facilitate its application in data sparse areas where only remotely sensed data can be obtained. Furthermore, in our definition, large scale does not limit the model to large sized rivers. This means the large scale model would be expected to include a substantial number of smaller streams a tributaries, with the minimum stream size determined by the application and data available rather than the model structure.
Approaches to regional or global scale river and floodplain simulation (Alfieri et al., 2013;Neal et al., 2012a;Paiva et al., 2011;Sayama et al., 2012;Winsemius et al., 2012;Yamazaki et al., 2011) often make a number of simplifications from the one-or two-dimensional shallow water models widely used at the reach scale (e.g. river lengths of 10-100's of km). Most of these simplifications are forced on the modeller due to insufficient information about the river geometry, channel roughness, floodplain topography and river discharge, and are compounded by computational resources that restrict the choice of model resolution, numerical complexity and ensemble size. River channel parameters such as depth, cross-section shape and friction are perhaps the most difficult to estimate on account of not being remotely observable (Alsdorf et al., 2007) and may therefore need to be estimated to obtain simulations of sufficient accuracy. Furthermore, current large scale inundation models are limited to using either no channels or simple rectangular ones, which is a situation that contrasts greatly with reach scale models where channel geometry can be complex and parameterised from survey data.
Despite a lack of data and suitable model structures, accurate simulation of channel flow and water level dynamics is essential to inundation simulation because even during large flood events the channel will still convey a significant proportion of the flow. Furthermore, a complex interaction between channel and floodplain flows might also be expected along the river plan form (Harwood and Brown, 1993;Trigg et al., 2012). As stated earlier, channel geometry data are unlikely to become available for many regions in the near future, meaning an alternative approach to using direct observations is needed. This is most likely to involve treating channel geometry as a model parameterisation, which can then be estimated along with friction. In this context, treating all channels as the same shape (i.e. rectangular) may introduce model structural errors that must then be compensated for, along with numerous other potential error sources, by calibrating the channel friction and depth parameters. By not allowing channel shape to vary, the friction and depth parameter will become more 'effective', which could lead to a number of spurious non-physical effects as friction and channel shape affect wave propagation in rather different ways and with different spatial signatures. At the reach scale, a number of studies have analysed the sensitivity of level and inundation simulations to model parameters, such as friction, and input data and have tried to identify the sources of uncertainty that are most prevalent (Apel et al., 2004;Domeneghetti et al., 2012;Pappenberger et al., 2005Pappenberger et al., , 2006. From these studies numerous modelling frameworks that account for uncertainty have been developed (Di Baldassarre and Montanari, 2009;Hall et al., 2011;Hostache et al., 2009). However, it is unclear to what extent the understanding of error sources and model sensitivity gained from these studies can be transferred to regional and global modelling, where the scale of the river network and observation data available differ. Previous large scale modelling studies have evaluated channel depth and friction as calibration parameters (Neal et al., 2012a;Paiva et al., 2013;Schumann et al., 2013), however the approaches lack the flexibility in model structure to incorporate more complex channel geometries. The method presented here will for the first time extend a large scale model to consider channel shape as a single continuous parameter (based on a power curve) using a unique method for determining channel wetted perimeter. This is the paper's fundamental novel conceptual advance, a model scheme that allows for structural changes to channel shape where computational efficiency is maintained while only introducing a single model parameter, thus maintaining a parsimonious model structure. Many existing models allow for complex channel geometries (e.g. Brunner (2010)) and trapezoidal channels, but the authors are unaware of any previous model that have taken this parameterised approach.
The overall aim of this work is to separate the effects of friction and geometry on wave propagation and river level such that more realistic large scale inundation models can be built and structural alternatives tested. However, there are a number of technical challenges to be overcome in order to be able to achieve this and this is what this paper describes. In subsequent sections of this paper, a hydraulic model structure is presented where the river channel depth, shape and friction can be described by three physically meaningful but continuous parameters, which might be calibrated or estimated from observations. Specifically, a method for approximating the channel shape using a power law is outlined, along with a novel and computationally efficient method for estimating the wetted perimeter of these channels. The approach to modelling the river channel is integrated within a two-dimensional inundation model (LISFLOOD-FP) using a simple-to-implement one-dimensional channel model (Neal et al., 2012a). After validating the method using a simple test case, the influence of channel shape friction and depth on river dynamics was investigated with the aid of a test case from a 60 km reach of the River Severn, UK.

Hydraulic modelling
To facilitate the development of the continuous channel geometry parameter a regional scale flood inundation model is needed. For this task, the hydraulic model of Neal et al. (2012a) and its approach to modelling sub-grid scale river channels within a two dimensional Cartesian grid floodplain inundation model was chosen. This model is implemented within the LISFLOOD-FP program because the authors are familiar with the code, however it should be possible transfer the proposed channel shape treatment to most one-dimensional hydraulic models because the variables being calculated are cross-section area and wetted perimeter. In the Neal et al. (2012a) model, river channels of any width below that of the floodplain model's cell resolution Dx and any bed elevation below that of the floodplain could be simulated. However, the channel geometry was assumed to be rectangular in order to keep the numerical scheme simple and computationally efficient. In this paper the model retains that functionality and computational efficiency, however an extension to the model is presented allowing channel cross-sectional geometry to be defined by a power function. This function relates the flow width w flow for a given depth of flow h flow to the bank-full width w full and bank-full depth of the channel h full using a single shape parameter s, such that the flow width for any flow depth below bank-full depth is defined as: The parameter s can take any value above 0 and produce a geometry. However, in physical terms, values below one will result in convex shaped banks, a value of one will lead to a triangular channel, and values above 1 will result in concave channels (e.g. a value of two is a parabolic channel, while the channel will tend towards trapezoidal then rectangular as s increases towards infinity). These geometries are illustrated by Fig. 1a. The channel is assumed to be symmetrical, which allows simple analytical relations to be derived between the flow width h flow and flow area A flow , given any value of s, as illustrated by Fig. 1b and Eq. (2): Bank-full width and depth are required as inputs to the numerical scheme, and are typically estimated from observations, hydraulic geometry theory or model calibration. The hydraulic radius (see Appendix A.1 for details) is defined by the flow area A t c;flow divided by the wetted perimeter P t c;flow , where flow denotes the flow area between two cells and c will be used from now onwards to denote a channel. For a rectangular channel the wetted perimeter is the channel width plus twice the minimum of the flow depth and the bank-full depth. However, for the power function shaped channel the wetted perimeter is determined by the function arc length. Unfortunately, calculating the arc length of a power curve is computationally expensive because there is no analytical solution, thus including this would not yield a computationally efficient model. In fact, even the analytical solution to the parabolic channel (s = 2) requires an expensive logarithmic function. Therefore, wetted perimeter was approximated using the flow width w t flow , which was calculated by Eq. (A5) and also shown in Fig. 1b, plus an additional component P h such that: The calculation and derivation of P h are described in Section 3. The numerical scheme is fully described in Appendix A and includes the equations for channel flow, floodplain flow, continuity and model time-step.

Derivation of wetted perimeter model
In Section 2 a simple method for approximating the wetted perimeter of an arbitrarily shaped channel was presented. This method requires the estimation of a wetted perimeter fraction P h , which can be added to the water surface widths to approximate the wetted perimeter of a channel of depth h and shape s. This section will present a linear regression method for estimating P h and demonstrate two key aspects of this approach: 1. That two quadratic models can be used to approximate P h to within 2% of that calculated using the trapezium method. 2. That the parameters of this model can be estimated for any value of s between 1.3 and 20 with a more complex regression model that need only be computed once at the beginning of a simulation.

Wetted perimeter approximation
The first step in approximating the wetted perimeter fraction P h was to derive a scaling such that an approximation could be derived for a channel of unit width. This involves scaling the values of depth of flow h, width of flow w and wetted perimeter fraction P h by the bank-full width w full , such that With this scaling w 0 will vary between zero for a dry channel and one for bank-full width or overbank flow. Scaled depth h 0 will vary between zero and infinity, but will typically be less than one because river channels are usually wider than they are deep. To account for deeply incised channels values of h 0 as high as two were considered. The second step of the process was to approximate the true wetted perimeter for the different shaped channels. To do this a near exact trapezium method was used to calculate the wetted perimeter for the channels with scaled depths of flow h 0 ranging for 0.01-2.0 sampled at intervals of 0.01, where channels have shapes from s = 1.3-20. This sampling accounts for any depth of flow between 0.01% and 200% of bank-full depth, channels more incised than this are assumed to be rare with accuracy likely to decrease gradually as h 0 increases beyond this range. Fig. 2a plots the scaled wetted perimeter P 0 for channels with shape s = 1.3, 1.5, 2, 5, 10, 20 against h 0 and illustrates how the increase in P 0 with depth is highly non-linear for small values of h 0 especially as s increases. Note that below s = 1.3 our method becomes inaccurate. Since this non-linearity would be very difficult to approximate with a polynomial, the scaled water surface width was taken away from the scaled wetted perimeter to create P 0;h , resulting in the relationship plotted in Fig. 2b (Note that water surface width is calculated as part of the model's momentum equation (Eq. (A5) in Appendix A), so it is already available at each time-step). The additional wetted perimeter component needed by Eq. (2) is therefore Fitting a polynomial to approximate P 0;h using linear regression was a trade-off between the order of the polynomial and the accuracy of the approximation, while partitioning the problem into smaller fractions of h 0 or using non fractional powers of h 0 could also increase the accuracy of the P 0;h estimates. Since the approximation needs to be implemented at every model time-step, there was a strong motivation for keeping the computational cost low. The use of fractional powers would make the wetted perimeter calculation almost as expensive as the momentum equation (Eq. (A2)), therefore the choice of polynomial was limited to integer power terms. Furthermore, P 0;h is still difficult to approximate with a single polynomial, therefore the data were split in three with h 0 between 0.05 and a threshold h 0TH modelled separately from values Testing different values of h 0TH indicated that a value of 0.2 provided a good compromise between P 0;h accuracy at high and low h 0 . The partitioned wetted perimeter fraction data are plotted in Fig. 2c and d. When apportioned in this way P 0;h could be approximated by a quadratic with a R 2 greater than 0.99 for all values of s between 1.3 and 20, as illustrated in Table 1. The approach could be made more accurate using additional polynomials, but this would add complexity and cost to the model for little benefit. The calculation of P 0h when 0:05 < h 0 6 h 0TH becomes where b values are parameters that depend only on s. When h 0TH < h 0 6 2 this equation becomes where P 0TH is the value of P 0h from Eq. (7) for threshold depth h 0TH .
For h 0 below 0.05 wetted perimeter is simply estimated given the top width (Eq. (A5)), which is a good approximation because the channels are concave in shape such that width increases rapidly  Table 1 Fit statistics for regression models of P 0;h and b parameters.
for small depths. Using high order polynomials (e.g. more b parameters) was investigated, however this was the simplest approach tested that could typically approximate P 0h to within 2% of the trapezium method estimates as discussed in the next section. The 2% accuracy is somewhat arbitrary, however the error is small relative to some of the other errors typically associated with hydraulic models, such as inflow discharge where rating curve errors up to 20% can be expected (Di Baldassarre and Montanari, 2009;McMillan et al., 2012). Using quadratics was also appealing because the function decreases to zero with h 0 , meaning the approximation of P 0;h would never be negative and always increase with h 0 because all of the b parameters are positive.

Approximation of b parameters
As stated previously the value of the b parameters depends only on s. Therefore, a model for how b varies with s is needed such that the b parameters can be calculated at the start of each simulation. Since this need only happen once at the start of the simulation the b parameter model can be significantly more complex than that used in Eqs. (7) and (8) where c parameters were found using a linear regression of each b parameter against s. Eq. (9) could be simplified and still provide good fits for some of the beta parameters, however since the computational cost of this model is negligible this was not pursued. The fit for each modelled b parameter against those found using the trapezium method is shown in Fig. 3 with R 2 values summarised in Table 1 parameter are presented in the supporting online material.

Accuracy of wetted perimeter approximation
To assess the accuracy of the channel approximation it was first necessary to understand the error introduced by assuming the wetted perimeter is equal to the water surface width (i.e. ignore P h ). For selected values of s, Fig. 4a plots the scaled depth against the percentage error in wetted perimeter when assuming the water surface width is the wetted perimeter. As all values of s are above 1 the channels initially widen faster than they deepen when h 0 is close to zero, while percentage error always increases with scaled depth. Errors are negative because the scaled wetted perimeter is always under predicted by the water surface width.
For selected values of s, Fig. 4b plots scaled depth against the percentage error in wetted perimeter estimates using Eqs. (7) and (8). The plot demonstrates that for s > 2 (e.g. channels with a parabolic to rectangular shape) the error in scaled wetted perimeter is less than 2% for all values of scaled depth. For values of s below 2 (e.g. channels between a parabolic and triangular shape) the percentage errors could be greater than 2% for scaled depth values below 0.05 with the greatest error in scaled wetted perimeter being 6.3% at scaled depth of 0.001 m and s value of 1.3. Errors in wetted perimeter for low s could be reduced by reducing the threshold scaled depth h 0TH at which the transition is made between Eqs. (8) and (7). However, this reduced the accuracy for higher values of s, which we decided were more important because most channels will have a more rectangular than triangular shape. The impact of channel shape approximation on wave simulation accuracy will be assessed using a simple test case in the next section.

Straight channel test
This section presents the results from a synthetic straight channel test case. The test has the three objectives: 1. Check for numerical stability when using the new wetted perimeter model. 2. Compare simulated levels from models with a range of s values (1.3, 2, 5 and 20) with the ''truth'' model for rectangular, triangular and parabolic channels, where analytical solutions exist for the wetted perimeter. 3. Conduct a one at a time sensitivity analysis of channel shape and Manning's n.

Test case description
The test case has a grid resolution of 500 m and a DEM for the floodplain and bank heights that slopes linearly from 100 m to 40 m above datum over a distance of 100 km. A discharge time series observed at the Saxons Lode gauging station on the River Severn (UK) between 1st June 2007 and 31st July 2007 was chosen and scaled by a factor of 2. Peak inflow for this event was 1625 m 3 s À1 and the original time series can be seen in Fig. 7. This hydrography was chosen because it includes both typical summer low flows and two flood events , which enabled a range of flow dynamics to be considered. This DEM and boundary condition setup was previously used by Neal et al. (2012). The model has a constant channel depth of 5.45 m, a width of 100 m and a constant Manning's floodplain friction of 0.035. Simulated levels, which are used to compare models with different parameterisations of n and s in the next section, were extracted at the centre point of the reach.

Accuracy of simulation using wetted perimeter approximation
In Section 3 the accuracy of the wetted perimeter approximation was established, however the numerical stability of the scheme, the impact of the approximation on depth simulation and the sensitivity of the model to channel shape relative to Manning's n will now be assessed. To assess the sensitivity of simulated depth to channel shape the analytical solutions for rectangular, parabolic and triangular channels were coded. The results from these models and the power curve shaped channel for four values of s (1.3, 2, 5 and 20) are plotted in Fig. 5a, while the difference between the approximate and analytical wetted perimeter for a parabolic shaped channel (s = 2) is plotted in Fig. 5b. As expected, these results show that depth simulated by power curve shaped channels falls between the triangular and rectangular models and covers the majority of the range of depths between the two. Mass balance errors were always at or below order 10 À8 for all simulations.
Regarding the sensitivity to s, maximum simulated depth varied from 7.24 m for s = 20-11.45 m for s = 1.3, a difference of 4.16 m. This would be a significant change in water level for a flood mapping application. To put this into context, Fig. 5c compares this range to the depths simulated by varying Manning's n from 0.015 to 0.065 in the analytical parabolic channel. This range of n values covers those found in almost all natural channels. The sensitivity of the model to n and s was similar for this test case. However, there are two key differences in the simulation dynamics when adjusting n or s. Firstly, changing n has a greater effect on the travel time of the flood waves for a given change in water level (a consequence of this is that the sensitivity to n might be expected to increase with distance). Secondly, low flow dynamics are relatively more sensitive to s than n, with high flow dynamics more sensitive to n for the channel parameterisations tested here. These results have established that the model is sensitive to both n and s, that these parameters affect model dynamics in different ways and that they will also have different characteristic length scales in space.

River Severn test
In this section the model is applied to a real world test case on the River Severn, UK. This site is subject to regular flooding, notably in 2000, 2007 and 2012, and has received extensive attention in the flood inundation literature focusing on both the flood hydraulics (Bates et al., 2006;Neal et al., 2011;Wright et al., 2008) and remote sensing of inundation (Garcia-Pintado et al., 2013;Giustarini et al., 2013;Mason et al., 2010;Matgen et al., 2011;Schumann et al., 2011). Using this site provides observation data for channel parameterisation and assessing the model  performance, while allowing an exploration of parameter sensitivity and model response.

Model setup
The hydraulic model was set up to simulate a 60 km reach of the River Severn between the Environment Agency gauging stations at Bewdley and Haw Bridge (Fig. 6). A 100 m resolution digital terrain model was created by aggregating 2 m resolution bare earth LiDAR digital elevation data from the Environment Agency (Fig. 6). The choice of resolution was motivated by computational cost alone and should not be seen as indicative of good practice. The embankment network between the gauging stations at Kempsey and Haw Bridge was subsequently added to the DTM elevations given location data from the UK National Flood and Coastal Defence Database and LiDAR elevation data (Fig. 6). Inflows to the model were based on gauge measurements from the seven Environment Agency gauging stations at locations shown in Fig. 6, with the downstream level set using observations of level from Ashleworth. The gauge discharges at the seven inflow stations are plotted in Fig. 7.
The locations of river channels were taken from UK Ordnance Survey mapping data, with the channel bank-full width, bed elevation (thalweg) and shape estimated at a total of 181 ground surveyed cross sections obtained from the Environment Agency. The greatest spacing between cross-sections was 500 m. To estimate channel shape a power function was fitted to the normalised relationship between channel width and depth for all cross-section measurements below bank-full depth. This is illustrated by Fig. 8 which shows the normalised width-depth relationships from the observation data and the fitted model for s. Residual errors demonstrated that for all cross-sections the error between observed and estimated depth was <10% of bank-full depth for >90% of the observations. Fig. 9 plots the river long-section for bed elevation, bank height, width and shape. Average channel shape in terms of s was 2.67 and between 2 and 5 for most of the reach. In a few locations the channel becomes more rectangular, with values close to 10. These values occurred mostly through the city of Worcester between 22 km and 28 km downstream or further upstream near a weir and lock gate, suggesting these are engineered sections. For simplicity the results section will use a uniform s value of 2.67 (mean shape), however this could be allowed to vary in space with width and depth and there is scope here for constraining the value of s in large scale models with geomorphic knowledge. Given that the typical cross-section spacing is greater than the model resolution, cross-section parameters were linearly interpolated to the 100 m model resolution between cross-sections. All other rivers in the domain (e.g. everything except the Severn main stem) were assumed rectangular and parameterised using hydraulic geometry in the same ways as described by Garcia-Pintado et al. (2013) and García-Pintado et al. (2015). Level simulation accuracy was not assessed on these rivers but we assume that they will accurately route flow to the Severn as demonstrated by Garcia-Pintado et al. (2013).
Finally it was necessary to include weirs in the model, the locations and parameterisation of these weirs were taken for an Environment Agency 1D model of the reach. Only the weirs at Tewkesbury and Diglis, both Crump weirs, are downstream of the water surface elevations used later for model evaluation and will influence levels at some upstream locations during low flow. The location of all weirs are included in Fig. 6, while details of the weir formulation used by the LISFLOOD-FP model can be found in the user manual . The model was run for a period of three months from 1st June 2007 to the 31st August 2007. The first 22 days of the 92 day simulation period were used as a model spin up period to ensure the model initial conditions, which were dry everywhere, had no influence on model performance during evaluation. Inflow data are plotted in Fig. 7.

Simulation sensitivity analysis
The model was not calibrated, however a series of simulations were run to establish the impact that changing Manning's n, bed elevation and shape of the Severn main stem has on the results. Channel Manning's values tested were from 0.015 to 0.07 at intervals of 0.005. Bed elevation z was varied globally by up to ±1.4 m at intervals of 0.2 m. Shape parameter s was not sampled systematically, instead the min, max, 10th, 90th, mean and median of the distribution of observed s was simulated. An additional simulation with s equal to 20 was also conducted to represent a rectangular channel. The simulation with a Manning's n of 0.03, mean observed s of 3.17 and observed bed elevation will be referred to as the first guess model. For this model, which has over 160 k cells, the computation time was 40.3 min on four cores of a 2.8 GHz Intel Xeon processor. An equivalent rectangular channel geometry model (note that this is not the s = 20 model) took 37.3 min for the same n and bed elevations. Time steps were as low as 5 s and typically below 10 s, indicating that the generalised channel geometry approach is computationally feasible relative to the rectangular channel model.
Model results were analysed at four gauging stations (Worcester, Diglis, Kempsey and Saxons Lode), with results at the three stations downstream of these (Mythe Bridge, Deerhurst and Haw Bridge) ignored due to their proximity to the downstream boundary at Ashleworth. Level simulations from the first guess model are plotted in Fig. 10 at the four analysis sites. The first guess model simulated both high and moderate water levels more accurately than the peaks and tended to under predict during recession from peak levels. These high flow errors are thought to result primarily from two sources (1) additional local runoff around Worcester and Kempsey not captured by a gauging station and (2) under estimation of peak flows at Bewdley due to rating uncertainty. This means a substantial component of this error is probably not due to the hydraulic model. Underestimation at Bewdley may have occurred because the flows at the Buildwas gauge 30 km upstream are up to 50 m 3 s À1 greater than at Bewdley. To investigate this issue the whole analysis presented in this section was repeated given Bewdley levels estimated based on a regression to historical levels at Worcester. Although this added nearly 200 m 3 s À1 to the flow to the river Severn and improved model performance, the conclusions from Figs. 11-13 remained the same. Therefore, results based on the published Environment Agency flow data are presented here.
To examine the sensitivity of simulated level to model parameters s, n and z, Fig. 11 plots the level response at Kempsey to each individual parameter. In all cases the solid dark line is the first guess simulation, with other lines representing the range of level simulations for that parameter. Here we can see that the model was most sensitive to Manning's n given the range of parameters tested and least sensitive to bed elevation. Sensitivity to n was expected, the similar sensitivity to s and the range of bed elevations tested (±1.4 m) also makes sense given the relative change in cross-sectional area (e.g. at bank full channel area is reduced by 1/3 for a parabolic channel relative to a rectangular channel, while the change in bed elevation of up to 2.8 m is approximately 1/3 of the mean observed channel depth). As parameter sensitivities are likely to be site, parameter range and performance metric dependent, the relative effect of each parameter on the model's dynamical behaviour is perhaps of more interest. Specifically, although all parameters can change levels, changing Manning's n had a greater effect on wave travel time and attenuation. This is illustrated by the results in Table 2 that describe the effect each parameter has on level at day 50 of the simulation period, during the rising limb of the largest flood event. The data also measure the range arrival time of the rising limb at Kempsey for the same event. The key result here is the relative sensitivity of level and timing to the parameters when normalised to channel shape, which shows that the Manning's n parameter is changing wave arrival time more than the other parameters relative to level. In fact, channel shape had the least effect on arrival time relative to level given this measure. On the falling limb Manning's coefficient appears to have even more of an impact on wave recession speed, however the behaviour is complicated by the presence of an increase in discharge around day 57.

Model performance
For each gauge the model performance, given 15-min resolution gauge data over the evaluation period, was calculated using the Nash-Sutcliffe efficiency measure. For this normalised statistic a value of 1 indicates a perfect fit between simulated and observed level and zero is equivalent to the mean level being used to predict the series. Below zero the score is no longer monotonic and indicates very poor model skill that would not be considered. Given the observed bed elevations the model performance at Kempsey given different channel shapes and friction is plotted in Fig. 12, which also includes an illustration of the channel shape formed by parameter s. From this we can make a number of observations regarding model performance. Firstly, whether calibrating channel shape, friction or both the maximum model efficiency is essentially the same given this measure. When the channel is parameterised with the mean of the observed channel shapes the optimal friction is somewhere between 0.03 and 0.035. If the channel shape becomes more triangular the optimal Manning's decreases, along with the maximum efficiency. This indicates that the loss of channel cross-sectional area for a given depth of flow can be partially compensated for by a reduction in friction. For a rectangular channel the optimal Manning's n is between 0.045 and 0.05, or 0.015 higher than for the channel with mean observed shape. The model performance is broadly similar and the parameter response surface has flattened relative to more triangular channel shapes (i.e. the sensitivity of model performance to n reduces as the channel becomes more rectangular). This would appear to indicate that the additional channel shape information is unnecessary, however the model with mean observed geometry calibrated to a friction far closer to the first guess for this reach. The increase in Manning's coefficient from 0.03 to 0.045 for the rectangular model will also delay the arrival of the rising limb around day 50 by $86 min over a distance of 30 km. From an analysis of all the observed events it was difficult to estimate the wave travel time from Bewdley due to the impact of tributaries, however this may be as low as 490 min, meaning the friction increase changed wave arrival time by up to $17% of travel time. Nash-Sutcliffe was also calculated over a number of time periods shorter than the full simulation duration. The lighter lines in Fig. 12 plot the model efficiency given only days 22-60 of the simulations period, which removes much of the low flow period towards the end of the simulation where discharge is not changing. In the absence of these low and relatively constant flows the response surface becomes sharper due to the loss of low levels that are predicted relatively well by all model parameterisations. The choice of model performance measure will therefore influence the model calibrations results, however the general conclusion of Manning's needing to be greater for the rectangular model was robust. Interestingly, even a poor guess at channel shape between the 10th and 90th percentiles of the observed shapes will move the optimum friction value closer to a first guess value. The other gauges show broadly similar results and are not reported here. Finally, for completeness, the model performance (Nash-Sutcliffe at Kempsey gauge) given bed elevation error and Manning's roughness is plotted for the mean observed channel and the near rectangular channel in Fig. 13. The results show the ability of both geometries to simulate water levels with similar accuracies, but that increased friction or bed elevation was required to obtain that optimum.

Discussion and conclusions
This paper has presented a simple and computationally efficient method for incorporating and testing structural changes in channel geometry into a large scale hydraulic model. The aim was to create a continuous variable that could be used to represent simple channel shape variability between near triangular to near rectangular. Analytical solutions were developed for the dynamic computation of channel top width and cross-sectional area, whereas wetted perimeters were approximated to maintain computational efficiency. Two analytical methods were tested, the first simplified the wetted perimeter to top width, while the second included an additional component to the perimeter based on two quadratic equations with parameters that depend on the channel shape. A simple test case was developed to test the robustness of the two approaches. For this test, the second method approximated wetted perimeter to within 2% of that estimated using a trapezoidal method for channels more rectangular than parabolic and always estimated wetted perimeter more accurately than using top width. It also had a negligible effect on simulation time, meaning it was adopted during subsequent tests. To our knowledge this is the first time such an approximation of wetted perimeter has been used by a large scale hydraulic model and therefore adds flexibility to the model structure that was not available previously. The approach provides a flexible modelling structure that can test different channel geometries as well as parameterisation schemes, and is a unique contribution for models that can run efficiently over large domains.
The new geometry model was validated using a simple test case based on a straight channel with constant channel properties and geometries that could be solved analytically, namely triangular, parabolic and rectangular channels. It was demonstrated that the approximation error for wetted perimeter had far less effect on simulated levels than the errors typically expected from real world applications of a hydraulic model. For a test case on the River Severn, the inclusion of additional channel shape detail, conditioned on surveyed cross-section data, was shown to improve model simulations of level given physically appropriate estimates of channel roughness and depth parameters. Calibrating channel Manning's in a rectangular channel could provide essentially the same performance in terms of Nash-Sutcliffe efficiency. However, this required a substantially higher Manning's coefficient (an increase of 0.015-0.02) above the more physically plausible value which negatively influence the simulated wave propagation time (by 1.4 h or $17% in this case). The insensitivity of the Nash-Sutcliffe efficiency measure to timing errors occurs because the discharge dynamics are slow in time relative to the timing errors introduced by calibrating to a higher Manning's coefficient. In a short reach like the one used here (30 km), the alteration of wave speed may be of limited consequence for practical applications at this scale, meaning that calibration of friction is still a useful tool for improving model performance as found by numerous studies (Horritt, 2000;Hostache et al., 2009;Mason et al., 2009;Pappenberger et al., 2005;Werner et al., 2005). However, as the reach lengthens, the cumulative effect of less (or non-) physically realistic friction in a calibrated rectangular channel simulation is likely to have an increasingly noticeable and negative effect on the accuracy of wave arrival time. Therefore, for large scale hydraulic modelling it may not be appropriate to simplify the channel geometry and then calibrate the model by adjusting the friction using an objective function based on level data (e.g. Neal et al. (2012)). On the other hand, large rivers that are wide relative to their depth may have sections that are effectively rectangular. In such cases the channel shape parameter may be near rectangular and be an unnecessary complication to the model structure that would ideally be avoided. Furthermore, for wide but shallow rivers with geometries closer to parabolic and triangular the wetted perimeter may not differ significantly from flow top width, meaning the scheme could be simplified in such cases. A simpler approach might assume top width is representative of wetted perimeter, but still calculate flow area because flow area will be different to the rectangular case by 1/2 for a triangular cross-section and 1/3 for a parabolic cross-section. The results also have implications for calibration schemes where friction is varied locally to fit water surface elevations (e.g. ) because, as a wave propagates downstream, a physically unrealistic increase in friction locally to fit water levels will need to be countered elsewhere if downstream wave arrival is to be correctly predicted. Specifically, these results demonstrate that shape can have a similar effect to friction locally, while having less influence on wave propagation than depth or friction changes, meaning it could provide an alternative model calibration parameter. The findings also showed that even a poor estimate of channel shape, defined as any shape between the 10th and 90th percentile of those observed, resulted in a calibrated channel Manning's coefficient closer to a physically plausible value that might be used in a reach scale study with observed cross-sections (i.e. less of an effective value). This is important because it suggests that including even basic knowledge of a channel's geomorphology in the model structure (e.g. Mersel et al. (2013)) could improve large scale hydraulic model calibration given observations of level, which are increasingly available from satellite altimeters on medium to large rivers (Maillard et al., 2015). Furthermore, in some cases it may be possible to use sequences of optical or radar imagery to get a first-order reach averaged estimate of the relationship between channel width and discharge (Smith and Pavelsky, 2008) above the lowest simulated discharge. Without the structural correction provided by the channel shape the model calibration may need to balance local water level simulation with longitudinal wave propagation and have difficulty simulating observed within channel level dynamics.

mDx ðA1Þ
where z c is the channel bed elevation, h t i is the cell water depth and Dx is the width of a cell. m is a meander coefficient which will be 1 for water surface slopes on the floodplain, but can be increased above 1 for rivers that meander within the grid and decreased below 1 for straighter rivers that are not aligned with the grid (e.g. the river length is shorter than the number of sub-grid channel cells multiplied by Dx). In this equation the water surface slope will be negative if the water surface in cell i + 1 is greater than i allowing flow to move in both directions and flow direction to change dynamically over a simulation. Given the water surface slope and SI units, flow between the two cells is calculated by where g is the acceleration due to gravity, Dt is the model time step, n is Manning's friction coefficient, subscript c indicates channel (subscript f will be used later to denote floodplain), and Q is the flow. The hydraulic radius R t c;flow and area of flow A t c;flow depend on the geometry of the channel and the depth of flow at the cell interface h t c;flow , which is defined as the difference between the maximum water surface elevation and bed elevation Given h t c;flow , the flow area and hydraulic radius are defined according to the equations in the next two sections.

A.1.1. Flow area equations
For a rectangular channel, multiplying the flow depth by the minimum channel width from the two cells (i and i = 1) is sufficient to obtain channel flow area A t c;flow . However for the power function shaped channel (Eq. (1)) the area of flow A t c;flow is defined as the minimum of the flow areas given the geometry of the sub-grid channel on each side of the cell interface such that because bank-full depth and area are known a priori. By taking the minimum area in this way the minimum flow area controls the hydraulics and therefore fluxes, while the possibility of having a large channel discharge to or from a cell with a smaller channel volume is avoided.
A.1.2. Momentum equation for the floodplain When the water depth in a model cell exceeds bank-full depth h i;full or the cell does not contain a sub-grid channel the model will simulate floodplain flow using the Cartesian grid and heights from a digital elevation model (DEM). Since these equations have been reported extensively by previous researchers de Almeida et al., 2012;Falter et al., 2014;Neal et al., 2011Neal et al., , 2012b only a brief explanation is provided here. To calculate floodplain flow the following equation is implemented at the cell interface where z f is the floodplain elevation and h full is zero where there is no sub-grid channel in a cell. For each cell interface the floodplain flow Q tþDt f ;iþ1=2 can be added to the channel flow Q tþDt c;iþ1=2 to calculate the total flow at the cell interface Q tþDt iþ1=2 , which will be used by the continuity equation to update cell volume over a time-step.

A.2. Continuity equation
The continuity equation is required to calculate the water depth h tþDt i;j in a cell given the change in cell volume from the previous time step, which results from flows calculated by the momentum equations (Eqs. (A2) and (A7)). Other source terms Q tþDt L;i;j that might include lateral inflows, evaporation, infiltration and rainfall can also be incorporated at this stage, such that, the updated cell volume is calculated using V tþDt i;j ¼ V t i;j þ Dt Q tþDt iÀ1=2;j À Q tþDt iþ1=2;j þ Q tþDt i;jÀ1=2 À Q tþDt i;jþ1=2 þ Q tþDt where V t i;j is cell volume and subscripts i and j reference the cell in two-dimensions. When the updated volume is less than bank-full volume V tþDt i;j < V i;j;full and the channel is rectangular, then the volume is converted to depth via Although presented in terms of volume, this is actually identical to the continuity equation presented by Neal et al. (2012a) except for the addition of the meander coefficient m. For the power curve shaped channel from Eq. (1) the volume is converted to depth by h tþDt i;j where the constant C i;j is calculated in advance using C i;j ¼ mDx i;j w i;j;full sðh i;j;full Þ À1=s s þ 1 For overbank flow when cell volume is greater than channel volume V t i;j > V i;j;full the equation simplifies to h tþDt i;j because bank-full depth is known a priori. The principal downside of the power curve shaped channel is the additional computational cost relative to the rectangular channel cause by the potentially fractional power functions in the equations for area (Eq. (A4)) and mass continuity (Eq. (A10)). Although this only affects within-bank flows since channel area and volume below bank-full can be pre computed.

A.3. Model time-step
As discussed by Bates et al. (2010) for an explicit scheme such as the one described above, the model time-step must be restricted to prevent numerical instability. The time-step constraint is defined by where max(H) is the maximum depth in the model domain and a is a coefficient below 1, which we set to 0.7. In this example the meander coefficient is the minimum of m and 1 (e.g. the meander coefficient only reduces the stable time-step if the river channel is shorter than the cell resolution).