Riverbed response to high flows with installed willow spiling using Multiquadric Radial Basis Function (RBFMQ)

To establish the success of any river engineering revetment, post-project monitoring is a critical part of the process. Vegetation structures often lack detailed biological or geomorphological data to aid management procedures and project planning. Apart from the survival and growth rates, data mapping detailed changes on the riverbed and its response to bankfull events are critical. Riverbed adjacent to two willow spilling structures on the River Stour in East Anglia, UK, was surveyed before and after high flow events. The random field data were gridded using Radial Basis Function with Multiquadric option (RBF MQ) to reconstruct the surface most accurately, as it was found to be the most effective out of three methods tested. Data were gridded, sliced and cross-sections were plotted and compared at each site before and soon after high flow events. Image maps were generated, and volumetric changes were calculated. At both sites, erosion dominated in the upstream and sedimentation in the downstream sections of the bed surveyed. Structure undercutting occurred at gravel site up to 29 cm which is a rate that needs a management intervention to prevent project failure.


Introduction
For observing small changes on river channel morphology, one of the most critical questions is how to, most accurately, represent the continuous 3D terrain model. Terrestrial and airborne laser scanning providing cloud datasets are now widely used in process based fluvial geomorphology (e.g. Pirotti et al. 2013, Clubb et al. 2017, Resop 2019), but can also show errors, especially where flowing water or vegetation occur (McKean 2009) and were shown as unsuitable for detailed river features (Wang et al. 2020). Thus, the remote data collection must be complemented with other field methods (Resop 2019) such as ground surveying. Here, a network of discrete, randomly spaced elevation data is manually collected. These are then reduced to a continuous function that generate a dense mesh of xyz data representing the model of the true topographic surface (Hardy 1971). The precision of any DTM depends on the accuracy of the approach used and conditions of the survey (Yilmaz 2007), but also on the data interpolation (gridding) function and its parameters (Watson 1999, Gao et al. 2020. Gridding is a function that interpolates field data to evenly spaced xyz grid nodes (Hardy 1971, Franke 1979) that can be shown as a surface plot and used for modelling and volumetric calculations. There is a range of gridding techniques available.
Radial Basis Functions (RBFs) are a family of interpolation methods and a primary tool to handle arbitrarily scarce data, to easily generalise to several space dimensions and to provide spectral accuracy in different types of applications (Wright 2003). RBF has been shown as the most effective group of methods to reconstruct smooth, manifold surfaces from point data and to repair incomplete meshes (Franke 1979, Carr et al. 2001. The method was used for most sophisticated CAD imagining in medicine, manufacturing, animation, architecture and other industries (e.g. Er et al. 2002, Michler 2011, Wang and Li 2011, Wee et al. 2011, ____________________ DOI: https://doi.org/10.33542/GC2021-2-01 Yoo 2011. RBF approximation techniques combined with machine learning algorithms such as neural networks are applied to optimize numerical algorithms (Skala 2011, Guang 2021 or for 3D surface reconstruction from scattered topography data (Hardy 1971, Carr et al. 2001, He et al. 2019, reconstruction of damaged images or inpainting removal (Liu et al. 2007). RBF can model complex surfaces from large data sets accurately and fast (Carr et al. 2001) and uses a single, automated, mathematical function applied to all the data points. The interpolation finds a complex shape that fits nearly exactly all the surveyed points and estimates elevation in places, where there is no data. Some of the other methods require either regularly spaced data or use existing data as estimates. The multiquadric RBF option (MQ) produced accurate models even with a small number of grid nodes (Ferreira 2003). This paper presents an example on how elevation coordinates from riverbed can be used and gridded to evaluate geomorphological changes, recorded over short periods of time, on small sections of the River Stour in East Anglia, UK, using RBF MQ. The results of this work are important for further management of the willow spiling structures installed at the river sites.

Study sites
Two soil bioengineering revetments using willow spiling to stabilise eroding riverbanks were built on the River Stour in East Anglia ( Fig. 1) and observed for changes on the riverbed first year after installation. This is known as the establishment period, which is critical for the project success (Allen and Leech 1997). The primary aim was to stabilise the eroding riverbanks by living structure instead of hard engineering (Anstead et al. 2012). Signs of significant loss of bed or bank material could point out gaps in the design that would gradually lead to project failure (Anstead and Boar 2010).
The River Stour is 108.5 km long and drains an area of 1,044 km 2 . Geology of the catchment is represented by Cretaceous chalk covered with boulder clay, fluvial sands and gravels. Soils along the river channel are represented by loam and clay floodplain soils naturally high on groundwater and freely draining slightly acid loamy soils. The land use is mainly arable land and horticulture with scattered settlements, grasslands and broadleaved woodlands. The highest discharges and prolonged wet periods occur during winter season. When the discharges are low, flows are supplemented by a water transfer scheme providing up to 3.8 m 3 /s of additional water into the river, more than double of natural river flow. The long-term annual rainfall in the watershed is 550 mm. Annually, around 40% of rainfall in the catchment enters the river and aquifer system and 60% is lost by evapotranspiration. The Base Flow Index (BFI), which describes the contribution of ground water to the river flow, varies between 52 and 43%.
The project sites were named according to the bank material as cohesive site C1 (sandy silt loam and loamy sand) and non-cohesive site N1 (sand and gravel). C1 is situated upstream from the N1 site, on 48.8 river km (from headwaters) and 24 m a.s.l., in the conservation area of water meadows near Sudbury, Suffolk. The discharge is natural plus the addition from the EOETS (Ely-Ouse to Essex Transfer Scheme), which brings water from the River Great Ouse from Cambridgeshire to supplement the low water supply in Essex (Anstead and Tovey 2017). The river channel at this site has been straightened in the past but is now naturally recovering its sinuosity. Meander radius is 29.3 m, amplitude is 127.3 m and site sinuosity is 1.15. Bankfull channel depth was 1.9 m and bankfull width 25.3 m. The site is also located on the main channel, about 100 m downstream of a weir. Banks are made of layers of silt with varying clay content and are typical with cantilever failures and slumping. Bed is silted and point bar is present at opposite riverbank.
N1 is situated 200 metres downstream from a confluence of flood alleviation and the old Stour channel near the village of Nayland in Essex. It is situated on a downstream end of a migrating meander with erosion up to 1.3 m/year. Meander radius is 52 m, amplitude 208 m and sinuosity 1.71. Banks are high but less steep than at C1. Bankfull depth is 2.3 m and width 22 m. Bed material is made of gravel and varies in form and depth. Both banks are also used for grazing and for rambling.
Grain size distribution was not performed for the surveyed riverbeds, but measurements from bankfoot were taken using a shear vane. C1 had significantly higher values of shear strength (44.8 ± 14 kPa) than N1 (7.30 ± 1.4 kPa).

Methods
Detailed survey of riverbed adjacent to spiling was performed before and after high flows that occurred in winter eight months after projects installation using TS Nikon DTM 330 (Fig. 2). Data were gridded using RBF MQ in Surfer 21 by Golden Software (2021). Out of 13 interpolation methods, Kriging and Natural Neighbour were also employed for comparison. Other methods were tried but they did not work as well for this dataset.
Kriging is a flexible method, quite effective with the default linear variogram for most data, although for larger datasets it can be slow. It is a precise or smoothing method and extrapolates grid values outside the data Z range (Golden Software 2021). It is a popular method as it produces appealing maps from irregular data (Yilmaz 2007). We used ordinary Kriging with linear variogram in our analysis.
Natural Neighbour is based on Thiessen polygons and can generate good grids from datasets in some but sparse data in other areas. It does not generate cell values in areas with no data. It is an exact interpolator and will not provide Z values beyond the Z input data range (Golden Software 2021).
Radial Basis Functions are a family of interpolation methods that are exact, which means that the interpolated surface will pass through each input data point (Hardy 1990, Yilmaz 2007. Each of the basis functions will produce a different interpolation surface. They work on a principle of smoothing parameter and the computations are based on the weighted combination of basic kernel functions to overlay the source data most accurately (Powell 1990). The kernel functions are analogic to variograms used in Kriging (Yilmaz 2007) and define the optimal set of weights to apply to the data points when interpolating a grid node (Golden Software 2021). The most used RBFs are Natural Cubic Spline, Thin Plate Spline, Gaussian, Multilogarithmic and Multiquadric methods (Hardy 1990, Yilmaz 2007, Balážovič 2015. RBF interpolation is a mesh-free method, meaning the nodes (points in the domain) need not lie on a structured grid, and does not require the formation of a mesh. It is often spectrally accurate and stable for the large number of nodes even in high dimensions. In a scatter point set, the calculation accuracy of multiquadric interpolation is strongly associated with the selection of shape factor, determined by the local point densities, therefore Adaptive Radial Basis Function (ARBF) interpolation algorithm was proposed (Gao et al. 2020). In our analysis, the shape factor R 2 has been calculated by default in Surfer as the squared length of diagonal of the data extent divided by 25 × number of data points (Yilmaz 2007).
The MQ method was presented by Hardy (1971) to solve a cartographic problem (Carlson and Foley 1991), to create a continuous function to reconstruct a topographic surface from elevation measurements (in Rocky Mountains National Park), (Wright 2003). This was found to be a function providing an exact fit of the data and a good approximation of the features of the surface (for example the position of peaks, depressions, valleys, junctions, or ridges). Additionally, it was automated and unbiased for generating contour maps of a topographic surface, which was a revolutionary approach for the time (Wright 2003). Although some analytical method of approximation existed, the resulting map was usually influenced by the individual cartographer. Furthermore, by having an automated function, mathematical methods could be used for the determination of landforms and their properties. The method has taken on new improvements and applications since (e.g. Carlson and Foley 1991, Carr et al. 2001, Goto et al. 2007, Sarra and Sturgill 2009, Deng et al. 2021. The original interpolation formula of the function was defined as (Hardy 1971): (1) where are the weights, x is the general n-dimensional vector and ( ) = (‖ − ‖) is the factor of distance (Skala 2011). Geometrically this compares to interpolating the data by a linear combination of n translates of the absolute value basis function | |, where the vertex of each basis function is located at one of the surface points (Wright 2003). However, this did not work well for peaks and valleys, therefore the absolute value basis function was replaced by a continuously differentiable one (Hardy 1972): where c is a determining constant (c ≠ 0), which will lead to a continuously distinctive interpolant. Instead of the absolute value basis function, this method compares to a linear combination of translates of the hyperbola basis function (Wright 2003). This method more accurately represented the real landform. Furthermore, instead of using a linear function of Euclidean distance basis function, Hardy (1972) proposed using a linear combination circular hyperboloid basis functions translated to be centred at each source point (Carlson andFoley 1991, Wright 2003): This function is infinitely differentiable ( ≠ 0), therefore techniques determining properties of the topographic surface can be applied. This was an excellent method for approximating a topographic surface from sparse, scattered measurements (Hardy 1986). It was not subjected to large oscillations such as the Fourier series methods and was able to account for rapid differences unlike the polynomial series methods (Wright 2003). Hardy (1986) named this new technique the "Multi-Quadric method" because he considered the principal feature of the method to be a "super-positioning of quadric surfaces". Further formulas were developed for 3-dimensional and multi-dimensional scattered data.
In addition to cartography, the method was effective in solving problems in hydrology, geodesy, photogrammetry, or geology. Franke (1979Franke ( , 1982 tested 29 available interpolation methods where RBF MQ provided best approximation in most tests (Wright 2003, Gao et al. 2020) and this method became more popular amongst mathematicians and researchers (Carlson and Foley 1991). RBFs were first introduced to solve problems of real multivariate interpolation (Frank 2014). As with other interpolation methods, RBF can be relatively inaccurate near boundaries and the boundary induced errors and correction techniques were addressed in the work by Wright (2003).

Fig. 2.
Scheme of the methodological steps from data collection in the field to establishing the changes using 2D cross-sectional or 3D volumetric data with RBF MQ. The output grid layer is a useful visual tool to determine the most critical locations.
The gridded XYZ data were cut to produce cross-sectional profiles through the generated gridlines. At each point, where the boundary line crossed a gridline, data point was generated. Cross section data were written to an ASCII data file (.DAT) that contained: (1) XY coordinates of the boundary line and grid line intersection, (2) Z value at the boundary line and gridline intersection, (3) accumulated horizontal distance along the boundary line and (4) boundary number used when more than one boundary line was contained in the file. The Y coordinate and the accumulated horizontal distance were then saved in a blanking file format (.BLN) and plotted as a cross section. Eight random riverbed cross sections were plotted through each surveyed area. They were positioned with the null distance close to the spilling and directed approximately perpendicular to the flow (Fig. 3). The plots are shown in the results section ( Fig. 4 and 5).
The interpolation with RBF is also a powerful tool for three-dimensional analysis. The cutoff volume between each of the grid surfaces and the horizontal plane given by the Z value for the highest full contour was calculated. The number of grid nodes used outside the blanked grid for the volume calculations was between 2,000 and 3,000. In Surfer, the volume is generated for each grid cell, and so the more grid cells available, the higher the accuracy of the final volume. The volume calculations were computed using three different methods: Extended trapezoidal rule, Extended Simpson's rule and Extended Simpson's 3/8 rule. If the results were close together, the true volume was close to those values (Hardy 1971). The net volume was then computed as the mean of the three values. To see, where along the elevation ranges the change happened, all elevation data have been used to compute percentiles and for colour relief plots. The whole process is summarised in Fig. 2.

Results and discussion
Riverbed adjacent to the two willow spiling structures was surveyed in November and in March using randomly spaced data points (Fig. 3). Eight control points were selected out of surveyed data to evaluate three gridding methods: Kriging, Natural Neighbour and RBF MQ. The Z values, the absolute differences between computed and the surveyed values, are shown in Tab. 2 and the best results are highlighted. Overall, Kriging provided the least accurate results. RBF MQ came as the most accurate with 5 out of 8 best results and NN achieved 4 best results. In the case of CP2 at N1 (Nov), both methods gave exactly same values, with only a 2 mm absolute difference from the original data. All methods performed better for control points in the middle of the surveyed field and worse near borders, which means that the selection of control points and the distribution of all data are important. Most significant error, 42 cm, was produced by Kriging in the case of CP1 at N1 (Nov), which came to almost 42 cm, but other methods also showed significant error: NN 32 cm and RBF MQ 21 cm which could mean some local irregularity on the riverbed. In summary for all control points, mean absolute difference for Kriging was 8.3 cm (13.8 cm), Natural Neighbour 6.2 cm (10.8 cm) and RBF with MQ 5.6 cm (6.7 cm). Overall, RBF MQ showed best results, however, more robust analysis with higher number of control points could be useful to confirm this conclusion. Results of previous studies testing gridding methods by Franke (1979), Ferreira (2003) or Yilmaz (2007) also showed RBF MQ as having the best outcomes. Accuracy of remote laser scanning is questionable in case of riverbed survey (Wang et al. 2020) and the resolution could be less than the changes observed. To report on the 2-dimensional analysis, 32 cross sections total were analysed. Notable changes were recorded at both field sites and individual plots are shown in Fig. 4 & 5. At C1, the morphological change within eight cross sections (CS) fell into the interval from -17% to +8%. All profiles recorded loss of bed material except two downstream profiles CS7 and CS8 (from -0.065 m 2 at CS3 up to +0.37 m 2 at CS8). At N1, changes were of a higher magnitude and fell into the interval -42% to +28% (from -0.351 m 2 at CS4 to +0.30 m 2 at CS8). Erosion prevailed on five upstream profiles and sedimentation occurred on remaining three downstream profiles CS6, 7 and 8, similarly to the previous situation at C1. Overall, computed from the cross-sections plotted, sedimentation prevailed at C1 (mean difference was +0.022 m 2  0.143 m 2 ) and erosion at N1 (-0.100 m 2  0.283 m 2 ).

Tab. 2. Areas (in m 2 ) of the random cross sections plotted through the gridded data files
before and after high flows at both sites. The percentage of difference in the cross-sectional area for the eight random cross sections at both sites was compared. Although there was significant erosion at N1 in the upstream section, the mean value has been influenced by the sedimentation in the downstream end. Hence Mann-Whitney test showed no statistical difference between the two sites (p = 0.442). Nonlinear correlation has been tested and the values for Pearson's product moment correlation computed. At both sites, there is a growing trend between the percentage of difference in the cross-sectional area and the distance along the spiling, emphasizing that the erosion prevails in the upstream and sedimentation in the downstream end of spiling as mentioned earlier. This could demonstrate effect the spiling may have on retarding river flow and sediment trapping. At N1, significant changes occurred also within the middle section of the spiling. To establish the immediate effect of the erosion processes on the spiling, the elevation changes at zero distance (closest to the spiling) were also compared. At the cohesive site, some erosion occurred, but it was not significant. The values varied between -1.7 cm and -0.93 cm and the mean value ±SD was -1.81 ±0.73 cm. At the N1 site, the toe scour was significantly higher, and values ranged between -29.14 cm to +3.27 cm. The mean change was -10.79 (±11.49) cm. Statistically, there was no significant difference between the two sites (p= 0.401) as again, the more extreme values at either end of the non-cohesive site influenced the position of the mean (Tab. 3).

C1
Cut-off plane (a plane taken through the highest common elevation) to determine the upper limit for cut-off (or fill-in) volume was chosen by the highest elevation present within the monitored area in before and after high flows (Tab. 3). Mean elevation over the C1 site has dropped by a small degree (by -1.3 cm) and slightly increased at the N1 site (by 0.78 cm). Cut-off volumes were calculated for this plane in Surfer and the mean value has been computed using the three rules described in the methods. Likewise, the mean elevation changes, there has been an increase in the cut-off volume at C1 indicating material loss of -0.227 m 3 , while there was an overall decrease in cut-off volume at the N1 site indicating 0.204 m 3 of material gain. This means that within the surveyed area of the riverbed, more material was deposited than it was eroded at N1.  To illustrate the changes within the elevation classes, the percentile distributions were analysed. Elevations have slightly decreased in all percentile groups at C1, with more extremes at values lower than 25 and higher than 75 percentile. At N1 the decrease occurred only within the 25 percentile of elevation and was very steep within the first 10% of samples. Elevations above 25 percentile increased slightly, meaning that the most obvious changes happened within the lower riverbed elevations. To illustrate this better, relief maps were plotted (Fig. 6) to indicate the spatial extend of erosion/accumulation processes. Higher elevations (in red) at C1 have slightly moved towards the lower elevations (in blue). At the downstream end, more space is occupied by the mid-range elevation. Although the slight colour shift shows slight erosion, this is not important in terms of undercutting or bank instability caused by bed scouring. At N1, the upstream section moved to lower elevation towards the bank and indicates area of bed scouring. On the other hand, deposition in the middle and downstream section is well illustrated by an increase in elevation. The accumulated material could have either originated from within the same site or from further upstream.
Considering the outcomes from the 3D analysis, over the surveyed area, more material was deposited at N1 site than it was eroded (+0.204 m 3 as opposed to -0.227 m 3 at the C1). On the other hand, results (both from 2D and 3D analyses) report on more extreme bed scour and undercutting at N1 site. The reason for this difference could be manifold, but important factor described earlier is the material which had very different engineering properties at each of the sites.
Non-cohesive, sandy material (N1) usually erodes more rapidly than cohesive one (C1) because cohesive particles such as clay minerals are chemically bound together, therefore are more resistant to flow (Hey 1975, Thorne and Tovey 1981, Thorne 1982, Lawler et al. 1997. The resistance of individual grains of sand and gravel to any movement is largely frictional, while clay particles have a high specific surface area and carry out electrical charges that will attract water molecules and cations providing for additional force known as cohesion (Craig 2004). The meander geometry of the N1 site, height of the bank and the previous riverbank retreat rates which were higher than on the C1 site could also contribute to this difference. What was common was, that in both instances, the bed erosion prevailed in the upstream section of the spiling and sedimentation occurred within the downstream sections. Although sedimentation in the downstream end may not have a significant adverse effect on the spiling, removing bed material in the upstream end at both sites can weaken the structure and, ultimately, can cause it to fail (Allen and Leach 1997). Ideally, spiling should extend between sedimentation zones and be sufficiently long to prevent any undercutting of weaker ends. The rate and timescale of scouring should be researched prior to installation because the bed scouring process may progress upstream. The removal of soil material from in front (and behind) the spiling has been listed as one of the most frequent reasons for project failure (Anstead and Boar 2010) and the spiling ends should be placed carefully to avoid areas with active bed scouring and high shear stresses (Allen andLeech 1997, Polster 2002).
Flooding and frequent bankfull events contribute to the bed material transport and by removing it can weaken the structure. During the observation period, November had exceptionally high monthly rainfall resulting in increased flows, February was a month with persistent high flows as the catchment was saturated. Peak high flow events occurred on 23. February (31.7 m 3 /s) and on 1. March (79.3 m 3 /s); values gauged at Lamarsh limnigraphic station, located on the river between the two research sites. The mean annual discharge for illustration was 2.68 m 3 /s. Spiling withstood well the flow events, but it was partially undercut. Furthermore, the sudden drop in water level weakened adjacent banks causing slumping upstream from the spiling. Recommendations on how to accommodate for high flows, bed and endpoint bank erosion when planning willow spiling project and other critical factors not dealt with in this paper, were reviewed by Anstead et al. (2012).

Conclusion
Getting the best model of the Earth surface is becoming less challenging with fast-developing quality of LiDAR technologies. In fluvial geomorphology, large-scale changes on river landscape and vegetation can be observed remotely and modelled with a sufficient accuracy. However, spatially small changes, especially those occurring under flowing water, are still challenging for laser scanning that needs to be complemented or replaced with traditional surveying methods. Here, the network of terrestrial data collected is still of good value if supplemented with an appropriate interpolation method.
The presented study demonstrated an application of Radial Basis Function with Multiquadric option for observing changes on riverbed before and after high flow season. This is not a common application as this method is more frequently used in other sectors such as medicine or manufacturing. A dense network of data was surveyed at two sites on the River Stour in East Anglia with installed vegetation bank stabilisation measures as part of the post-project monitoring process. Randomly collected data were processed and interpolated using Kriging, Natural Neighbour and RBF MQ. The last method showed the best results in a test which also confirmed the conclusion from previous studies. RBF MQ was used to create a dense mesh of data which allowed for subsequent 2D and 3D analysis. Colour relief maps were produced to identify the most critical parts that could lead to potential project failure and the volumetric analysis reflected the amount of material eroded.
This pointed out exactly the scale and location of retreat or sedimentation processes on the riverbed. In both cases, erosion dominated in the upstream and sedimentation in the downstream sections. Particle size distribution of riverbed was not carried out but shear strengths in the bank foot were measured, giving considerably higher values at C1 than at N1 site. Bank and meander geometry also played a part but with only two sites it was not possible to make any conclusions.