High-magnitude flood events are among the world’s most widespread and significant natural hazards and play a key role in shaping river channel–floodplain morphology and riparian ecology. Development of conceptual and quantitative models for the response of bedrock-influenced dryland rivers to such floods is of growing scientific and practical importance, but in many instances, modeling efforts are hampered by a paucity of relevant field data. Here, we combined extensive aerial and field data with hydraulic modeling to document erosion, deposition, and vegetation changes that have occurred during two successive, cyclone-driven, extreme floods along a 50-km-long reach of the bedrock-influenced Sabie River in the Kruger National Park, eastern South Africa. Aerial light detection and ranging (LiDAR) data and photography obtained after extreme floods in 2000 and 2012 (discharges >4000 m3 s–1) were used to generate digital elevation models (DEMs) and provide the boundary conditions for hydraulic modeling (flow shear stresses for three discharges up to 5000 m3 s–1). For the Sabie River study reach as a whole, DEM differencing revealed that the 2012 floods resulted in net erosion of ∼1,219,000 m3 (∼53 mm m–2). At the subreach scale, however, more complex spatial patterns of erosion, deposition, and vegetation change occurred, as largely controlled by differences in channel type (e.g., degree of bedrock and alluvial exposure) and changing hydraulic conditions (shear stresses widely >1000 N m–2 across the river around peak flow). The impact of flood sequencing and relative flood magnitude is also evident; in some subreaches, remnant islands and vegetation that survived the 2000 floods were removed during the smaller 2012 floods owing to their wider exposure to flow. These findings were synthesized to refine and extend a conceptual model of bedrock-influenced dryland river response that incorporates flood sequencing, channel type, and sediment supply influences. In particular, with some climate change projections indicating the potential for future increases in the frequency of cyclone-generated extreme floods in eastern southern Africa, the Sabie and other Kruger National Park rivers may experience additional sediment stripping and vegetation removal. Over time, such rivers may transition to a more bedrock-dominated state, with significant implications for ecological structure and function and associated ecosystem services. These findings contribute to an improved analysis of the Kruger National Park rivers in particular, but also to growing appreciation of the global diversity of dryland rivers and the relative and synergistic impacts of extreme floods.
In an era characterized by rapid environmental change and variability, increasing research attention is being directed to the role of extreme hydroclimatic events such as storm rainfall, flooding, and drought in the shaping of Earth’s surface. High-magnitude floods are among the world’s most significant natural hazards, and they play a key role in the shaping of riparian environments across a wide range of physiographic and hydroclimatic zones (Woodward et al., 2010). Drylands (hyperarid, arid, semiarid, and dry subhumid regions) are one of the most extensive hydroclimatic zones, covering 41% of Earth’s surface and sustaining 36% of the world’s population (United Nations, 2016). Drylands are characterized by net annual moisture deficits resulting from low annual precipitation and high potential evaporation, and typically by strong climatic variability. Although precipitation regimes vary widely, many drylands are subject to extended dry periods and occasional intense rainfall events. Consequently, dryland rivers are commonly defined by long periods with very low or no flow, interspersed with infrequent, short-lived, larger flows. On any given river, this results in a high ratio of large to small flows (McDermott and Pilgrim, 1983), resulting in highly skewed flood frequency distributions, and regional and relative flood frequency curves that are usually steep, as the slopes are typically established by a few, very large floods (Tooth, 2000). Some of these very large floods may be of sufficient magnitude or impact to be termed “extreme” (e.g., Gupta, 2000) or “catastrophic” (e.g., Thompson and Croke, 2013).
For many dryland rivers, highly variable flow regimes are key to their morphological development (e.g., Tooth and Nanson, 2011; Tooth, 2013) as well as the maintenance of important riparian habitats (e.g., van Coller et al., 2000; Kingsford, 2006; Parsons et al., 2006; Stromberg et al., 2007; Sandercock et al., 2007; Jaeger et al., 2017). In some physiographic settings, variable flow regimes and diverse riparian vegetation assemblages combine with limited sediment supplies and heterogeneous bedrock morphologies to produce dryland river morphologies and dynamics that differ markedly from fully alluvial rivers, particularly those in humid temperate regions (van Niekerk et al., 1995; Heritage et al., 1999, 2001; Wohl and Achyuthan, 2002; Tooth and McCarthy, 2004; Jansen, 2006). In recent decades, greater research focus has been directed toward such “bedrock-influenced,” “bedrock-controlled,” or “mixed bedrock-alluvial” dryland rivers (Heritage et al., 1999, 2001; Tooth et al., 2002, 2013; Tooth and McCarthy, 2004; Keen-Zebert et al., 2013), but morphological, sedimentological, hydraulic, and ecological data remain limited. This paucity of data hampers efforts to develop conceptual and quantitative models of the morphological, sedimentary, and vegetative response of these types of dryland rivers to past, present, and future climatic changes, including the importance of shifts in flood frequency-magnitude relationships, flood timing, and flood sequencing. These data and knowledge gaps are becoming increasingly significant as drylands are now widely considered to be some of the regions most vulnerable to future hydroclimatic changes (Obasi, 2005; IPCC, 2007; Wang et al., 2012), and this limited model development restricts attempts to develop environmentally sound, sustainable management practices for such rivers.
Some previous work, however, provides a basis for improved model development. For instance, morphological responses in some bedrock-influenced dryland rivers have been shown to be controlled by rare, high-magnitude floods that are responsible for doing the most geomorphic work (Baker, 1977), contrasting with many temperate rivers, where the bankfull flood with a 1–2 yr return period has typically been viewed as being the dominant channel-forming discharge (Wolman and Miller, 1960). Furthermore, cross-section–scale models of quasi-cyclical channel development have been proposed both for subhumid Australian (Nanson, 1986) and semiarid South African (Rountree et al., 2000) rivers, whereby channels undergo net deposition during lower flow periods, but extensive erosion occurs during high-magnitude floods that leads to system “resetting.” To date, there has been little extensive testing of these conceptual models using field or modeling data, and the possible influence of climate-driven changes to flood regimes remains poorly known. In addition, there have been few detailed considerations of the variation in channel responses over longer river reaches that may include a variety of channel types or starting states.
Nevertheless, generating the requisite field data to improve models of bedrock-influenced dryland river response is not a trivial exercise. Many of these rivers are located in sparsely populated areas, where instrumental flood records are either absent or of limited length, gauges often are destroyed during high-magnitude flood events, and scale and/or morphological variability can introduce considerable spatial complexity. Nonetheless, significant recent advances have been made to our understanding of bedrock-influenced dryland river response to past floods through the use of sediment archives (Zawada, 1994, 2000; Macklin et al., 2010; Woodward et al., 2010). In addition, sophisticated remote survey technologies such as light detection and ranging (LiDAR) and increased computer-processing capabilities have opened up the possibilities of capturing high-resolution topographic data and embedding these data in more sophisticated morphodynamic models to characterize river responses to recent and potential future floods (Milan and Heritage, 2012; Croke et al., 2013; Thompson and Croke, 2013; Baggs Sargood et al., 2015). In recent years, computational modeling has led to significant insights into the flow and sediment dynamics and consequent morphological responses of fully alluvial rivers (Nicholas, 2010; Nicholas et al., 2013). To date, there have not been similar advances in the understanding of bedrock-influenced river dynamics, where sediment supply limitations, resistance to scour, and complex roughness and flow partitioning can have significant influences on erosion, deposition, and resultant channel, bar, island, and floodplain development.
Against this backdrop, this paper demonstrates how a combination of field investigations and improved data capture and processing capabilities can be used to quantify the geomorphic impacts of high-magnitude floods on bedrock-influenced dryland rivers. We focus on the Sabie River, one of several large (>100 km long), bedrock-influenced rivers in the Kruger National Park in eastern South Africa (Fig. 1), where recent cyclones have generated extreme floods (defined here as flows >4000 m3 s–1). The aims of the paper are to: (1) combine aerial and field data with hydraulic modeling to characterize and explain erosion, deposition, and associated vegetation changes that occurred during cyclone-driven extreme floods in January 2012; (2) compare the changes that occurred during the 2012 floods with the impacts of the January/February 2000 extreme floods; and (3) synthesize these findings to refine and extend a conceptual model of bedrock-influenced dryland river response, including accounting for the role of event sequencing in determining channel and vegetation dynamics. Beyond the Sabie and other Kruger National Park rivers, our findings have more generic relevance for improving understanding of the dynamics of bedrock-influenced dryland rivers, both in South Africa and farther afield.
The Sabie River drains a 7096 km2 catchment that straddles the border between Mpumalanga Province, South Africa, and southern Mozambique (Fig. 1). The headwater tributaries begin in the Drakensberg to the west (∼1600 m above sea level [asl]), and the river descends rapidly onto the lower-relief Lowveld (∼400 m asl) and Lebombo zones (∼200 m asl) in the east (Fig. 1B). The middle reach lies within the boundaries of the Kruger National Park (Fig. 1C). Annual average rainfall is highest in the uplands (∼2000 mm yr–1) and declines rapidly toward the South Africa–Mozambique border (450 mm yr–1), where it is exceeded by average annual potential evapotranspiration (1700 mm yr–1). Rainfall occurs mainly in the austral summer (November through March) and normally results from convective thunderstorms, although occasional, high-intensity rainfall events can result from cyclones that form over the Indian Ocean and track inland. For example, the maximum daily rainfall recorded for Skukuza (Fig. 1C) between 1912 and 2001 was 103.5 mm (Kruger et al., 2002), while 2 to 5 day rainfall totals can exceed 200 mm (Heritage et al., 2001). The flow regime of the Sabie River reflects this hydrological regime, with summer floods being separated by long periods of winter low baseflows (generally <50 m3 s–1) that are supplied from dolomitic aquifers in the western headwaters. Although water abstractions have altered the low-flow regime, flood flows are unaffected, and the river remains unimpacted by engineering structures or other human activities over considerable lengths within the national park. As such, the Sabie River represents an excellent example of a near-pristine river that can be investigated to improve our understanding of bedrock-influenced dryland river and vegetation dynamics.
The Sabie River is underlain by a variety of sedimentary, metamorphic, and intrusive and extrusive igneous rocks (Fig. 1B). Long-term incision of the Sabie River into these heterogeneous lithologies has generated a wide “macrochannel,” within which narrower channels, bars, islands, and floodplains occur (Fig. 2). Lithological differences influence the longitudinal profile of the Sabie River (Fig. 1B) and thus affect the hydraulic conditions, the distribution and thickness of alluvial sediment, and channel morphology (Cheshire, 1994; van Niekerk et al., 1995). In the middle reach, extensive low-relief plains border the incised macrochannel, and sediment supply is mainly from tributaries and within-channel sources (e.g., alluvial bars and islands) rather than from hillslopes. Sediment is dominantly sand and fine gravel (median grain size 1–2 mm), as derived from weathering and erosion of the gneisses and leucogranites that underlie significant portions of the river (Fig. 1B; Heritage and van Niekerk, 1995; Broadhurst et al., 1997). Repeated marked changes in channel morphology occur as the distribution and thickness of sediment over bedrock change.
Heritage et al. (1999, 2004) defined and mapped six principal channel types along the middle Sabie River, leading to the definition of subreaches A to V based upon the dominant channel type present (Figs. 1D and 2; Table DR11): mixed single thread, mixed braided, uncohesive mixed anastomosed, cohesive mixed anastomosed, bedrock anastomosed, and mixed pool-rapid. In all cases, the “mixed” descriptor refers to a “mixed bedrock-alluvial” state and indicates the presence and influence of bedrock in an otherwise alluvial channel type (Heritage et al., 1999). The distinction between “uncohesive” and “cohesive” reflects the degree of consolidation of sedimentary deposits, generally distinguishing between deposits composed of comparatively mobile sand and fine gravel, and deposits characterized by a significant component of silt and clay, which tends to impart a greater resistance to erosion (Heritage et al., 1999). Each of these channel types is associated with various morphological units, including different bar forms (e.g., lateral, point, lee), islands, and floodplains (Fig. 2).
On the Sabie River, channel types and their diverse morphological units are important influences on vegetation patterns (van Coller et al., 2000). Elevation in relation to the macrochannel is associated with vegetation pattern variations, indicating that flood frequency and water availability are key determinants. For instance, species such as Ficus sycomorus are found along the floor of the macrochannel, while Spirostachys africana is typically found higher up on the macrochannel bank, reflecting changing proximity to the water table and flood frequency (Birkhead et al., 1996). Nevertheless, variations in soil, substratum, and nutrient status, and differences in bedrock exposure and morphological unit distribution complicate this simple elevation-vegetation association (e.g., van Coller et al., 2000). In particular, the degree of alluviation is also a key factor, with species such as Combretum erythrophyllum being closely associated with the mixed braided channel type, and Phragmites mauritianus associated with mixed pool-rapid sites. In contrast, Breonadia salacina and Syzgium guineense are more closely associated with bedrock anastomosing sites, where bedrock-core bars form the fundamental morphological units. Bedrock exposure provides opportunities for anchoring, allowing the establishment of these species, which are unable to root firmly in alluvium (van Coller, 1993). In addition, the microscopic seeds of B. salacina become trapped in cracks in the bedrock, increasing the potential for germination (van Coller et al., 2000).
CYCLONE DANDO FLOOD EVENT
In mid-January 2012, Cyclone Dando impacted on southern Africa, leading to widespread heavy rainfall (450–500 mm in 48 h). Many rivers draining into and through the Kruger National Park experienced extreme floods, resulting in widespread erosion and sedimentation (see Heritage et al., 2015). These floods occurred just 12 yr after Cyclone Eline also resulted in extreme rainfall and flooding in Kruger National Park (January/February 2000). For instance, between 5 and 10 February 2000, a total of 544 mm of rainfall was recorded at Graskop in the upper catchment, and a total of 245 mm was recorded at Skukuza in the middle catchment (Heritage et al., 2001). Along a 101-km-long reach of the middle Sabie River, peak discharge during the resulting 2000 floods ranged between ∼4000 and 7000 m3 s–1 (Heritage et al., 2001), and the associated erosion and deposition led to the transformation of channel types in many subreaches, including changes both to more alluvial and less alluvial states (Heritage et al., 2004). Widespread removal of woody riparian vegetation also occurred (Pettit and Naiman, 2006). Between the 2000 and 2012 floods, only low to moderate flows (<1000 m3 s–1) occurred (Fig. 3), and significant erosion and deposition were limited. This situation provides an opportunity to examine the relative and synergistic impacts of two extreme floods occurring in sequence.
We developed rating equations for eight sites (Fig. 1D, T1 to T8), based upon the relationship between site-specific simulated water levels and discharges. These rating relations were then used to estimate site-specific discharges based upon real-time kinematic (RTK) global positioning system (GPS) measured strand-line data. Discharge estimates for the Sabie River ranged between 4470 m3 s–1 and 5630 m3 s–1, with the peak estimates located ∼8 km downstream of the Sand-Sabie River confluence (Figs. 1C and 3). Overall, and in line with other lines of physical and anecdotal evidence (see Methods section), these modeling results suggest that the 2012 floods did not exceed the peak stage or extent of the 2000 Cyclone Eline floods, which reached ∼6000–7000 m3 s–1 toward the downstream end of the study reach (Heritage et al., 2001, 2004).
This study employed a combination of aerial image acquisition and analysis, digital elevation model (DEM) generation, field survey, and hydraulic modeling. Aspects of our methods and analytical approaches are similar to those employed in other studies of landscape response to extreme hydroclimatic events, including storm rainfall and flooding (e.g., Baggs Sargood et al., 2015; Tseng et al., 2015).
Following the Cyclone Dando floods, aerial LiDAR data and photography were obtained on 30 May 2012 for 50 km reaches of three rivers (Sabie, Letaba, and Olifants) in the Kruger National Park. Southern Mapping Geospatial surveyed the three rivers using an Opetch Orion 206 LiDAR, and a Rollei Aerial Industrial Camera (AIC) with a 60 megapixel P65+ Phase One digital Charge coupled device (CCD), flown from a Cessna 206 at 1100 m altitude. This paper focuses on the imagery obtained from the Sabie River study reach (Fig. 1C; Milan et al., 2018a), where aerial LiDAR data and photography were also available from 2004, having been collected several years after the large floods of January/February 2000. Given that no major floods occurred between 2004 and 2012 (Fig. 3), the 2004 imagery and LiDAR serve as a benchmark for assessment of the subsequent impacts of the 2012 floods. Following the approach in a previous study (Heritage et al., 2004), the 50 km study reach was subdivided into subreaches A to V, reflecting the dominant channel type present prior to the 2012 floods (Fig. 1D; Table DR1 [see footnote 1]).
DEM Generation and Error
From the 2004 and 2012 Sabie River LiDAR surveys point cloud data for ground returns were used to produce DEMs. The 2004 LiDAR point density was 45,869 points km–2, and in 2012, average point density was 409,318 points km–2. Aerial photographs taken coincident with the LiDAR survey also allow for a visual assessment of morphological and vegetation cover change. A DEM of difference (DoD) that accounted for error in each of the two surveys was then produced to allowed assessment of the spatial variations in the magnitude of erosion (sediment loss) and deposition (sediment gain) in each subreach and the 50-km-long study reach as a whole. The 2004 DEM was also used to provide the boundary conditions for two-dimensional (2-D) hydraulic modeling of the 2012 floods. In addition, point elevation data for vegetation returns were used to produce a DoD for the vegetation, enabling assessment of flood-related changes to vegetation cover, with changes confirmed through a visual inspection of aerial imagery captured before and after the floods.
Field Surveys of High-Stage Indicators
To provide an input to the hydraulic modeling (see following), in May 2012, RTK GPS surveys of high-stage indicators were undertaken along select ∼300–500-m-long sites along the middle Sabie River within the Kruger National Park (Fig. 1D, T1 to T8). Choice of survey sites was dictated by access and safety considerations, particularly the presence of dangerous animals. Despite the 4 months that had elapsed between the January floods and the surveys (a time gap imposed by the availability of funding), strand lines of organic debris (e.g., branches, twigs, reeds) were generally well preserved. This finding is in line with other studies that have shown how trim lines, scour marks, and flotsam can be well preserved in dryland environments and so can be used to reconstruct flood hydraulics and hydrograph characteristics, even some considerable time after the flood (e.g., House and Pearthree, 1995; Greenbaum et al., 1998). In our study reach, the fresh condition of the debris and identifiable “best before” dates on some of the embedded plastic bottles showed clearly that these strand lines were from the January 2012 floods. We acknowledge that the receding flood can leave several strand lines, depending on local conditions, but because surveys and subsequent data review focused on capturing the elevations of the highest strand lines at a given site, these surveys provide an indication of the highest stage reached by the 2012 floods. At the Low Level Bridge crossing near Skukuza (Fig. 1C), a roadside marker indicates the limit of the 2000 floods, and it stands at a higher elevation than the strand lines from the 2012 floods. At this location, therefore, the 2012 floods were not as large as the 2000 floods. Anecdotal accounts from park rangers suggest that this finding applies more widely along the middle Sabie reach, and it is supported by the absence of any damage during the 2012 floods to the tarred road running adjacent to the macrochannel margins along the right bank, as this road had been extensively damaged in the 2000 floods.
2-D Hydraulic Modeling
As noted already, the 2004 LiDAR data for the Sabie River were used to provide the boundary conditions for hydraulic modeling of the 2012 floods. Horritt and Bates (2002) noted that many of the roughness factors represented by the roughness coefficient in one-dimensional (1-D) models are integrated into the modeling process in 2-D models. As such, a nominal Manning’s n roughness value of 0.03 was used in JFLOW, a 2-D depth-averaged flow model. JFLOW is a commercial 2-D flow modeling tool noted for its ability to handle large data sets through the use of a graphics processing unit–based computation. JFLOW was developed as a solution to harness the full detail of available topographic data sets such as those available from LiDAR, and to investigate overland flow paths (Bradbrook, 2006). Simplified forms of the full 2-D hydrodynamic equations are used in the model, but the main controls on flood routing for shallow, topographically driven flow are captured (Bradbrook, 2006).
We undertook a series of steady-state simulations for the Sabie River as nine connected subreaches across a range of flows (350, 3500, and 5000 m3 s–1). The 350 m3 s–1 flow is the approximate magnitude of many annual floods (Fig. 3), while the higher flows are more characteristic of the periodic large or extreme floods. The aim was to find a simulated discharge that produced water-surface elevations close to those measured from the RTK GPS–surveyed strand lines in the field, so providing peak flow estimates for different subreaches along the middle Sabie River. The spatial patterns of shear stress at different discharges also provided information on the potential for channel morphodynamic responses and riparian vegetation impacts.
To validate the modeling, comparisons were made between the simulated water-surface elevations and the RTK GPS–surveyed high-stage indicators for sites along the Sabie River, examples of which are shown in Figure 4A. Simulated water-surface elevations generally matched the strand lines very well (typically ±0.25 m), giving confidence in the hydraulic performance of the model. Modeled and surveyed flood inundation extents were also generally well matched (Fig. 4B). The highest simulated flow of 5000 m3 s–1 is slightly higher in elevation at site T3 compared with the strand-line measurements (Fig. 4A), suggesting that JFLOW slightly overestimates discharges at this location. Farther downstream at sites T5 and T6, however, simulated water surfaces plot lower than measured strand lines (Fig. 4A), suggesting slight underestimation of discharges, and possibly indicating the importance of inflows from the Sand River upstream (Fig. 1C). This suggests that during the 2012 floods, peak discharges were slightly in excess of this simulated 5000 m3 s–1 flow.
Comparison of the available aerial imagery for the 50-km-long study reach of the Sabie River suggests subdued morphologic development in the period 2000 through to early 2012. As such, the 2004–2012 DoD represents change caused primarily by the January 2012 floods. For each of the subreaches A to V along the middle Sabie River, the DoD (Fig. 5) illustrates the patterns of erosion and deposition, and the sediment budget (volume/area) quantifies the magnitude of erosion and deposition (Fig. 6). An indication of the hydraulic conditions during near-peak flows is provided by the modeled shear stress distributions for the 5000 m3 s–1 flow (Fig. 5).
Downstream Patterns of Fluvial Geomorphological Change
The DoD, based on the difference between the DEMs (Fig. 5), reveals ∼3,344,000 m3 of erosion and ∼2,125,000 m3 of deposition over the 50 km study reach, indicating net erosion of ∼1,219,000 m3 (53 mm m–2) during the 2012 floods. The sediment budget (Fig. 6) reveals broad patterns in erosional losses and depositional gains, with erosion more dominant in the upstream ∼25 km reach and deposition more dominant in the downstream ∼25 km reach. Patterns and magnitudes of erosion and deposition, however, are not consistent within or between the dominant channel types in each subreach (Fig. 6), and there is no simple or consistent correlation with near-peak shear stress distributions (Fig. 5). Subreaches A (mixed single thread) and B (mixed pool-rapid) showed an approximate balance between erosion and deposition (Figs. 5 and 6). In these subreaches, there were alternating zones of high shear stress (up to 900 N m–2) interspersed with zones of lower shear stress (<500 N m–2; Fig. 5). Subreaches C to K contained a wide range of channel types and showed mainly net erosion, although the short subreach E (bedrock anastomosed) was an exception, as this had an approximate balance between erosion and deposition (Figs. 5 and 6). The most notable zones of erosion were evident in subreach F (cohesive mixed anastomosed), especially in the 3 km immediately downstream of the Sand River confluence, the downstream half of subreach I (mixed braided), and the upstream third of subreach K (uncohesive mixed anastomosed; Figs. 5 and 6). In all these subreaches, shear stresses showed considerable spatial variation (Fig. 5), with no obvious spatial correlation with local morphologic change. Subreaches L to Q all showed net deposition, regardless of channel type and near-peak-flow shear stress distributions (Figs. 5 and 6). In subreaches R to V, net erosion was again dominant, regardless of channel type or near-peak-flow shear stresses, with the exception of subreach T (mixed pool-rapid), where there was net deposition (Figs. 5 and 6).
Closer examination of subreaches D to K in the vicinity of the Sand River confluence revealed finer-scale erosional and depositional patterns, which enabled a more detailed investigation of the relationship to changing hydraulic conditions. This 15-km-long sequence of subreaches is of particular interest because it includes five of the six channel types present along the middle Sabie River (Figs. 1D and 2; Table DR1 [see footnote 1]). Steady-state 2-D hydraulic simulations are presented for the three different discharges of 350, 3500, and 5000 m3 s–1 (Fig. 7). At 350 m3 s–1, all subreaches have extensive zones with low shear stress values (typically <250 N m–2). Nonetheless, some subreaches show considerable spatial variability; for instance, in subreach F (cohesive mixed anastomosed), shear stresses are locally high (up to 900 N m–2) around a prominent island in the upstream part of the subreach, relatively low (<200 N m–2) on the outside of the subsequent bend in the macrochannel, and higher (>300 N m–2) on the inside of the bend (Fig. 7A). As expected, the higher flow simulation (3500 m3 s–1) results in a general increase in shear stresses, but especially in subreaches D, E, F, G, and H, which include zones of relatively high shear stress of 800–1100 N m–2 (Fig. 7B). The highest flow simulation (5000 m3 s–1) reveals a further general rise in shear stresses, but it generally maintains the broad spatial patterns of relatively low and high shear stress zones evident in the 3500 m3 s–1 simulation (Fig. 7C). More widespread zones of high shear stresses (800–1200 N m–2) are evident in subreaches D and E, and in the anabranches on either side of large alluvial islands in subreach F (Fig. 7C). The high shear stresses in subreach E are not associated with any significant erosion (Figs. 5 and 6), possibly in part because of the limited alluvial cover in this bedrock anastomosed channel type, but erosion does show a clear association with peak shear stresses (∼1100 N m–2) in subreach F, and this is particularly noticeable in the two anabranches that diverge around the alluvial island at the downstream end of this subreach (Figs. 5, 6, and 7C). Shear stresses up to ∼1600 N m–2 are found toward the downstream end of subreach G (bedrock anastomosed) and the upstream end of subreach H (cohesive mixed anastomosed), and these reaches are characterized by net erosion (Figs. 5, 6, and 7C).
Local Riparian Vegetation Change
The channel morphological changes and shear stress distributions (Figs. 5–7) help to explain patterns of vegetation loss through the study reach. The DEM of vegetation difference for the full study reach is shown in Figure 5. Notable vegetation losses occurred at the upstream end of subreach I and the downstream end of subreach K, in subreach N, in subreach Q, and at the upstream end of subreach U (Fig. 5). Subreaches E, J, and V did not experience any significant vegetation loss (Fig. 5). Local patterns of vegetation removal and survival were highly complex, however, as illustrated by an examination of the changes in the vicinity of the Sand River confluence (Fig. 8). The DEM of vegetation difference for subreaches D to K indicates that during the 2012 floods, partial stripping of late successional vegetation occurred (Fig. 8E). For example, a swath of mature riverine woodland was removed from the right-hand side of the large alluvial island at the downstream end of subreach F (cohesive mixed anastomosing) and the upstream end of subreach G (bedrock anastomosing). Such areas of vegetation loss appear to be associated with zones of high shear stress attained during the higher flows (Fig. 8C). In general, the exposed edges of islands appear to be more susceptible to the loss of large trees, compared with the better-protected island interiors. Very few areas that had well-established vegetation were stripped to bedrock.
Local Hydraulic, Topographic, and Vegetative Influences on Erosion and Deposition
Figures 5 through 8 show that clear, consistent, cause-and-effect relationships among shear stress distributions, erosional and depositional magnitudes, and vegetation changes are not apparent, suggesting that the observed flood-related morphological changes are likely an outcome of multiple factors. Along the Sabie River, erosional and depositional dynamics are undoubtedly complex; for instance, areas of alluvial cover could be stripped to some depth on the rising limb of a flood but covered with sediment on the falling limb, possibly resulting in net deposition over the course of the event. Furthermore, erosion magnitudes are also limited in many locations by the presence of bedrock outcrop or outcrop that lies at shallow depth beneath thin alluvial cover.
To investigate in greater detail the potential influences of local hydraulics, topographic roughness, and vegetation upon the magnitudes of erosion and deposition, we focused on a 500-m-long section of the dominantly cohesive mixed anastomosing subreach F downstream of the Sand River confluence. Using the difference between the shear stress grids for the 5000 m3 s–1 and 350 m3 s–1 flow simulations, we determined at-a-point magnitudes of shear stress change representative of: (1) shear stress increase, indicative of that observed on the rising limb of a flood and representative of energy gain and erosion potential; and (2) shear stress decrease, such as observed on the falling limb of a flood and representative of energy loss and depositional potential. These data (∼77,000 points) were then plotted against corresponding local erosion (negative values) and deposition (positive values) for each 2 m grid pixel on the DoD (Figs. 9A and 9B). Figure 9A shows that there is no straightforward relationship between shear stress increases and the magnitude of erosion. Erosion of up to several meters of sediment thickness occurred even in areas that showed very small shear stress increases (Fig. 9A), probably representing erosion of unconsolidated sand. Peak erosional losses, in the region of 4–6 m thickness, occurred where shear stresses locally increased by ∼300–700 N m–2 (Fig. 9A), for at these shear stresses, even relatively cohesive deposits may be susceptible to erosion. Notably, the extent and magnitude of losses tended to decrease after this peak (Fig. 9A), indicating that the areas with the greatest erosion were not necessarily associated with the zones of highest local shear stress increases. Visual inspection of the aerial imagery suggests that the zones that experienced local shear stress increases of between 2000 and 5000 N m–2 were locations with only limited initial alluvial deposits and greater bedrock exposure, with erosional losses thereby moderated by local sediment availability rather than available erosive energy. Figure 9B shows that there is also no straightforward relationship between shear stress decreases and the magnitude of deposition. The greatest deposition (>2 m thickness) is evident where shear stresses locally decreased by 200–400 N m–2 (Fig. 9B). The number of point locations with larger shear stress decreases is less, and this accounts for the apparent reduction in the magnitude of deposition evident toward the right-hand side of the graph (Fig. 9B).
Figure 9C shows the relationship between detrended elevation (representative of local topographic roughness) and erosion and deposition magnitudes. Again, no straightforward relationship is evident; up to several meters of erosion or deposition occurred across almost the full range of topographic variability (–4 to +4 m; Fig. 9C). Nevertheless, the topographic highs tended to have a much wider spread of erosion and deposition magnitudes, with high points tending to show a dominance of erosion (Fig. 9C). The plot seems to suggest a lower envelope curve, whereby areas of the bed with greater protrusion (increasingly positive values) were subject to greater erosion, possibly because they were areas with greater exposure to flow or had the thickest preflood sediments.
The presence of vegetation is known to influence fluvial processes, and especially patterns of erosion and deposition (Bywater-Reyes et al., 2017), so Figure 9D shows pre- and postflood LiDAR vegetation height returns for the 2004 and 2012 DEMs plotted against erosion and deposition magnitude. The hypothesis here is that the locations where vegetation removal has taken place may show an association with erosion, while locations where vegetation has survived may show an association with deposition. No clear relationship between vegetation heights and erosion and deposition magnitudes is evident, although there is a slight tendency for the areas with no or low vegetation height (particularly <2 m) to have been associated with greater erosion and deposition (Fig. 9D). More erosion (typically up to 5 m, possibly down to bedrock) is also evident compared to deposition (typically up to 2 m) for areas with vegetation heights <2 m. This may indicate that this vegetation was more susceptible to stripping during this event. Lower, less-scattered erosion and deposition magnitudes tended to occur in areas where more established (taller) vegetation was present, for example, on some of the alluvial islands in subreach F. This tendency is seen most clearly for the post-2012 flood vegetation height returns. Higher erosion and deposition magnitudes are evident in areas that had taller vegetation height in the 2004 imagery, possibly reflecting greater sediment mobilization in areas experiencing vegetation removal during the 2012 floods.
The results from the aerial image acquisition and analysis, DEM generation, field survey, and hydraulic modeling (Figs. 5–9) provide an opportunity to analyze and interpret flood-related change along a bedrock-influenced dryland river in greater detail than has been possible previously. The data sets also provide scope for comparison with the impacts of the January/February 2000 extreme floods on the Sabie River.
Changes during the 2012 Floods
During the 2012 floods, morphological change did not occur with any consistency within or between the dominant channel types in each subreach (Fig. 6), and there was no simple or consistent correlation with near-peak shear stress distributions (Fig. 5). In most subreaches, the magnitudes and patterns of erosion and deposition during the 2012 floods were not sufficient to transform the dominant channel type. Only three of the 22 subreaches showed visual evidence of a change in channel morphology, namely, subreaches C, I, and P, all of which showed a switch from dominantly mixed braided to dominantly mixed single thread (Fig. DR2 [see footnote 1]). For subreach C and I, this switch was associated with net sediment loss, but subreach P showed a net sediment gain (Fig. 6), likely owing to the development of lateral bars and floodplains (Fig. DR2 [see footnote 1]).
Close inspection of the aerial photographs for the wider study reach shows that some subreaches had more complex assemblages of channel types than is indicated by the simple subreach division according to dominant channel type (Fig. 1D; Table DR1 [see footnote 1]), with both cross-stream and downstream variations in morphologic complexity evident. For instance, subreach F (Fig. 8A) is characterized by a mixed assemblage of morphologic units associated with both uncohesive and cohesive mixed anastomosing channel types. Bedrock control becomes more evident toward the very downstream end of this subreach, with bedrock anastomosing then becoming the dominant channel type in subreach G (Fig. 8A). The greater the bedrock exposure within a subreach, the greater is the tendency for diversity in the assemblage of channel types, and thus within any given subreach, erosional and depositional patterns will reflect the patchy losses and gains from the various morphological units (e.g., lateral bars, point bars, lee bars, islands) that comprise the channel types present (Fig. 2). These patchy losses and gains in large part will be related to local hydraulic changes during floods, especially areas of shear stress increase and decrease (Figs. 9A and 9B).
Comparison of the 2000 and 2012 Floods
The nonlinear and spatially variable patterns of erosion, deposition, and morphological change in the 2012 floods are generally consistent with the response reported for the larger 2000 floods, with all channel types being subject to varying degrees of alteration in the distribution and thickness of sediment over bedrock (Rountree et al., 2000, 2001; Heritage et al., 2004). In the 2000 floods, however, morphological change tended to be more consistent across the preflood channel types. Some channel types (e.g., bedrock anastomosing) remained essentially unchanged, but there was more widespread transformation of other channel types to more alluvial or less alluvial states. In the 2000 floods, aerial imagery suggests that net erosion dominated over the middle Sabie reach as whole, with many cohesive mixed anastomosed subreaches in particular experiencing alluvial “stripping” to become bedrock anastomosed. In the 2000 floods, no single hydraulic parameter demonstrated a strong and consistent correlation with channel change (Heritage et al. 2004), and just as in the 2012 floods, change appeared to have been controlled by a complex combination of factors, including spatial and temporal variations in flow energy levels and associated sediment transport, and tributary inputs of water and sediment. During the 2000 floods, for instance, the Sand River is estimated to have delivered in excess of 65,000 tonnes of sediment, leading to increased sedimentation in the lower part of the study reach (Heritage et al., 2004). A similar pattern is evident following the 2012 event, with net deposition occurring in subreaches L through Q (Fig. 6), most probably reflecting sediment input from the Sand River. Aerial photograph comparisons show that most of this sediment was derived from reworking of alluvial bars and islands in the middle reaches of the Sand River.
This comparison of the 2000 and 2012 floods also provides insights into the role of flood sequencing and relative flood magnitude in determining vegetation changes along the Sabie River. For instance, although the 2000 floods were able to strip some of the subreaches hosting cohesive mixed anastomosed channel types, fragments of late-stage successional (“mature”) vegetation survived on remnant islands, as well as on some other older alluvial units. In these and other subreaches, however, some sediment stripping and removal of the mature vegetation occurred in the subsequent 2012 floods (e.g., subreaches F through I; Fig. 8E), even though these floods were smaller. An interval of 12 yr dominated by low to moderate floods (Fig. 3) appears to have been insufficient time for recovery of cohesive mixed anastomosed and other channel types to their pre-2000 flood condition, and additional partial stripping occurred because remnant sediments and surviving vegetation were left exposed to the full force of the subsequent floods, rather than being in the protective interior of larger, more coherent alluvial landforms (Heritage et al., 2015).
Rivers that have undergone extreme flood disturbance events provide opportunities for documenting and analyzing how channel-floodplain morphologies and associated riparian vegetation assemblages develop over space and time. There have been many studies of the impacts of historic high-magnitude floods in dryland rivers or wet-dry subtropical rivers (see review in Tooth, 2013), including comparisons of the differential impacts of floods on closely adjacent reaches or rivers (Thompson and Croke, 2013), but far fewer studies of the impacts of sequences of historical floods (although see Huckleberry, 1994; Fryirs et al., 2015). Against this backdrop, the investigations along the Sabie River provide a rare opportunity to examine in detail the relative and synergistic impacts of two extreme floods occurring in sequence, and to develop models of bedrock-influenced dryland river development.
Flood-Related Erosional and Depositional Patterns
Few studies have investigated in detail the relationships between flow hydraulic parameters and flood erosion and deposition in large, bedrock-influenced dryland rivers. For the Sabie River 2000 extreme floods, Heritage et al. (2004) investigated the effects of flood slope, shear stress, and stream power upon morphological response but failed to support the correlations reported in earlier research undertaken on rivers in different physiographic and hydroclimatic contexts (e.g., Howard and Dolan, 1981; Wohl, 1992; Wohl et al., 1994; Benito, 1997). No single hydraulic parameter demonstrated a strong, consistent correlation with erosional and depositional patterns on the Sabie River at the reach or subreach scale (Heritage et al., 2004). Focusing on shear stress (Figs. 5, 7, and 8), the findings of this study support these broad conclusions, while highlighting how a combination of high-resolution topographic data and hydraulic modeling nonetheless can provide insights into the localized, patchy patterns of erosion and deposition that occur during such extreme floods (e.g., Fig. 9).
Flood-Related Vegetation Dynamics
Previous studies of dryland rivers, both within the Kruger National Park and farther afield, have shown the importance of riparian vegetation for restricting erosion during large or extreme floods and/or enabling postflood “recovery” (e.g., Baker, 1977; Osterkamp and Costa, 1987; Sandercock et al., 2007; Pettit and Naiman, 2006; Tooth, 2013). In many studies, emphasis has been placed on the role of vegetation in facilitating channel, island, and floodplain development during the relatively low-magnitude floods that occur during “building” phases (e.g., Schumm and Lichty, 1963; Burkham, 1972; Osterkamp and Costa, 1987; Lisle, 1989; Hooke and Mant, 2000; Rountree et al., 2001; Greenbaum and Bergman, 2006). In the southwestern United States, groundwater–surface water interactions, influenced by catchment management activities, have also been shown to facilitate riparian vegetation growth (Webb and Leake, 2006). While the conditions and process thresholds (e.g., velocities, shear stress) that give rise to riparian vegetation removal during floods in dryland rivers still warrant further investigation (Thornes, 1994), flood-related vegetation losses have been studied in relation to catastrophic floods in small limestone streams in central Texas (Baker, 1977), as part of the “arroyo cycle” phenomena in alluvial streams in the southwestern United States (e.g., Graf, 1983; Hereford, 1984, 1993), and in sand-bed streams in the Great Plains of the United States (Friedman et al., 1996a, 1996b).
The comparison of the 2000 and 2012 Sabie River floods provides valuable additional insights, particularly by highlighting how partial—and potentially complete—loss of vegetation may occur in phases linked to the temporal clustering of large or extreme floods. The evidence from subreaches F to I (Figs. 8D and 8E), for instance, shows that vegetation that survived one large or extreme flood may be vulnerable to removal during a successive flood that occurs within a couple of decades, even if that flood is of lower magnitude. This can be attributed to the limited time available for the sediment deposition and woody vegetation development that is a necessary precursor for island and floodplain rebuilding, and to increased exposure of remnant sediments and vegetation to the full force of subsequent floods.
Models of Bedrock-Influenced Dryland River Development
In common with studies of some other dryland rivers (discussed earlier herein), early investigations in the Kruger National Park (e.g., Heritage et al., 2001; Rountree et al., 2001) outlined how the rivers are characterized by alternating phases of building and stripping. During long periods of quiescent low-magnitude floods, sediments build up within the bedrock macrochannel. Vegetated bars develop at the macrochannel margins, or islands develop between multiple anabranches, gradually reducing macrochannel cross-sectional area. The vegetative root network plays an important role in strengthening these bars and islands, essentially by lowering the potential for alluvial stripping. In addition, mature trees are able to capture large quantities of organic material and propagules from upstream, and this is likely to help perpetuate successional development of riparian vegetation at the morphological unit scale (cf. Gurnell et al., 2001; Pettit and Naiman, 2006; Merritt et al., 2010). During rarer, higher-magnitude floods, however, shear stresses can exceed the resistance thresholds of the accumulated sediments and vegetation, resulting in stripping back to the bedrock macrochannel template. Although broadly correct, these early analyses tended to consider mainly changes at the “whole-system” scale, rather than the subreach changes to individual channel types. In addition, although volumes of erosion and deposition and time scales of building/stripping remained largely unquantified, there was an implicit assumption that complete (or near-complete) stripping occurs during the high-magnitude floods.
By providing insights into the relative and synergistic impacts of two successive extreme flood events, including the localized patterns and volumes of erosion and deposition at the subreach scale, this study builds on this early work and the subsequent studies by Heritage et al. (2004, 2015). Here, we propose a further refinement to previous conceptual models of bedrock-influenced dryland river development in the Kruger National Park by incorporating flood sequencing, channel type, and sediment supply influences (Fig. 10). This model incorporates the initial starting state and the relative magnitude of successive floods, and it reflects the fact that following periods of sediment accumulation, various degrees of erosion can take place during high-magnitude floods, including partial stripping that leaves remnant sediment and vegetation on the bedrock template. Specifically, the temporal development from a relatively low-gradient, “mature” alluvial channel starting state (Fig. 10A), and a higher-gradient, “young” bedrock starting state (Fig. 10B) is dependent upon the relative magnitude of flood events (low, moderate, large/extreme) and the frequency and order of these events.
For both starting states, an increasing frequency of low-magnitude floods will promote sediment deposition and vegetation colonization (Fig. 10). For the lower-gradient alluvial starting state, occasional moderate floods may not initially cause significant morphological changes (Fig. 10A), although some uncohesive deposits may convert to cohesive deposits as silt is draped over the surface and vegetation becomes more established. An increasing frequency of moderate floods, however, may initiate partial stripping of sediment and vegetation (Fig. 10A).
For the higher-gradient bedrock starting state, an increasing frequency of low-magnitude floods promotes sediment deposition, but at a slower rate in comparison to the alluvial starting state, owing to higher energy levels and more limited vegetation cover (Fig. 10B). An increasing frequency of moderate floods may also allow some uncohesive sediment to accumulate, some of which may become more cohesive over time as silt drapes develop and vegetation establishes (Fig. 10B). In these moderate floods, flow energy is not capable of fully stripping sediment from the macrochannel, enabling vegetation to mature gradually, and sediment that is stripped may also be replaced by new sediment supplied from upstream.
In both types of channels, however, an increasing frequency of large or extreme floods will initiate stripping in alluvial channels (Fig. 10A) or maintain the channel in a largely bedrock state (Fig. 10B). For the alluvial starting state, widespread partial stripping of cohesive deposits may take place across the macrochannel, but if sediment supply is high, uncohesive sands may be deposited on the falling limb on these partially stripped surfaces. Such within-flood patterns of erosion and deposition occurred in parts of the Sabie River study reach during the 2012 floods (see also Knight and Evans, 2017). Following a partial stripping event, subsequent floods promote channel development along one of two different pathways (i.e., redeposition or further erosion), depending upon factors including the size of successive floods, the length of time between these floods, and additional factors affecting vegetation recovery (e.g., successional competition, herbivory damage). Relatively long gaps between large or extreme floods (e.g., 20–30 yr), with intervening smaller floods and favorable conditions for vegetation re-establishment, may enable recovery or transformation to a more alluvial state, particularly in subreaches that have been only partially stripped. Where the time gap is relatively short, however, such as between the 2000 and 2012 extreme floods, recovery or transformation to a more alluvial state is unlikely, and the trend is likely more toward net stripping down to a bedrock template. In subreaches F through I on the Sabie River, for instance, cohesive mixed anastomosed channel types were partially stripped during the 2000 floods, and remnant sediments and vegetation were left exposed and experienced further losses during the 2012 floods (Fig. 8). We propose that once stripping has been initiated, and other large or extreme floods occur in close succession, then the trend is commonly toward further stripping. Although complete stripping did not occur along the Sabie River during the 2012 floods, it did occur farther north on the Olifants River (Milan et al., 2018b). With an increasing frequency of large or extreme floods, complete stripping may eventually occur along the Sabie River; with removal of increasing volumes of the sediment stored in bars and islands, initially alluvial subreaches will be transformed to largely sediment-free bedrock subreaches (Fig. 10A).
The historical evidence from the Sabie River and other Kruger National Park rivers demonstrates that it does not take more than one or two large or extreme, closely spaced floods to remove many meters of alluvium that had accumulated within the bedrock macrochannel. Furthermore, optically stimulated luminescence ages show that most alluvium within the macrochannels is no older than a few hundred years (Heritage et al., 2015), which indicates that alluvial stripping occurs regularly in these extreme flood-prone systems, and a recovery period on the order of a few hundred years is required to attain a fully alluvial state with late successional riparian forest. Given that some climate change projections highlight the potential for a southerly shift in cyclone tracks and increased landfall over South Africa and Madagascar (Fitchett and Grab, 2014) and an increase in rainfall quantities during wet seasons (MacFadyen et al., 2018), it is possible that the Sabie River will no longer experience the prolonged periods of lower-magnitude floods that are necessary to enable the “end-member” alluvial state to develop. Instead, the Sabie River and other Kruger National Park rivers may experience a state change to a more bedrock-dominated system, with near-complete stripping of cohesive sediment occurring, and only thin veneers of unconsolidated sediment being distributed across the macrochannel. This will have significant impacts on the ecological structure and function of these rivers. In particular, recovery of some vegetation species will be restricted by the limited availability of alluvial substrate and increased exposure to the impacts of clusters of large or extreme floods. Essentially, along these supply-limited, flood-prone river systems, alluvial sediments can be conceptualized as being only in temporary storage, and with every stripping episode, incremental erosion of the underlying bedrock template will occur. Bedrock lowering rates are unquantified, but the wide vertical and horizontal joint spacing (typically >0.5 m) in the gneisses and leucogranites that underlie much of the study reach (Fig. 1B) precludes extensive hydraulic plucking, and so in any individual high-magnitude flood, bedrock lowering is likely to be negligible (i.e., millimeter-scale). Nevertheless, by this process, the rivers will continue to etch their macrochannels into the landscape, thereby contributing to overall landscape denudation.
Using a combination of aerial imagery acquisition and analysis, DEM generation, field surveys, and hydraulic modeling, this paper characterized and explained channel and vegetation response to two successive, cyclone-driven, extreme floods along a bedrock-influenced dryland river system in the Kruger National Park, eastern South Africa. During the 2012 floods on the Sabie River, net erosion of ∼1,219,000 m3 occurred over the 50-km-long study reach, although individual subreaches experienced varying degrees of erosion and deposition. By contrast with the 2000 extreme floods, there was only limited evidence of channel type switching, despite significant sediment reworking and redistribution along the study reach. Locally, partial stripping of mature vegetation occurred, with the margins of islands left exposed by the 2000 floods being particularly susceptible to further erosion during the 2012 floods.
By synthesizing the findings regarding flood-related erosion/deposition patterns, channel change, and vegetation dynamics, we have presented a new conceptual model for bedrock-influenced dryland river development. This model builds upon previous conceptual models of river development in the Kruger National Park by outlining how different pathways of channel development and associated vegetation dynamics depend upon the initial channel state (alluvial or bedrock), changes in flood magnitude-frequency relationships, and flood sequencing. Although developed for the Sabie River and other Kruger National Park rivers, many of the general concepts may be transferable to other rivers where periodic extreme flood events and riparian vegetation interactions are known to be key geomorphic drivers, including other large, bedrock-influenced dryland rivers, both within South Africa and farther afield. Examples include rivers in arid and semiarid regions of the United States, such as central Texas (Baker, 1977), as well as rivers in subtropical regions subject to alternating extremes of above- and below-average rainfall and flooding, such as southeast Queensland, Australia (Croke et al., 2013).
In recent years, considerable advances have been made in computational modeling of alluvial river dynamics, but similar advances in the understanding of extreme flood-impacted, bedrock-influenced rivers has been slower. Use of a similar combination of data sets and methods as outlined in this paper and other recent studies (e.g., Thompson and Croke, 2013; Baggs Sargood et al., 2015) will support model development and thus help to address this gap in knowledge. For drylands in particular, which are widely considered to be some of the regions most vulnerable to future hydroclimatic changes (Obasi, 2005; IPCC, 2007; Wang et al., 2012), this may have significant implications for improved assessment of potential changes to the geomorphology and ecology of riparian zones and associated ecosystem services.
This project was funded through the Natural Environment Research Council Urgency Grant NE/K001132/1. We would like to thank SANParks for supporting this research. We are grateful to SANParks for providing the 2004 LiDAR and aerial imagery data. The raw 2012 LiDAR and aerial imagery data may be accessed from Milan et al. (2018a). JBA Consulting is acknowledged for allowing use of JFLOW. The paper was improved as a result of reviews by Vic Baker and three anonymous reviewers, and the comments of Bruce Shyu (Associate Editor) and Brad Singer (Science Editor). We would like to thank the University of Hull for funding Gold Open Access for this paper.