A Computer Aided Approach for River Styles—Inspired Characterization of Large Basins: A Structured Procedure and Support Tools

Abstract

This paper presents a systematic procedure for developing a characterization and classification of river reaches inspired by the River Styles Framework, through which insight can be gained about the understanding of river behavior. Our procedure takes advantage of several computer based “tools”, i.e., algorithms implemented in software packages of various types, from “simple” Excel sheets to sophisticated algorithms in Python language, in general all supported by Geographic Information Systems (GIS). The main potentially useful, existing tools for this specific aim are discussed here, revealing their strengths and weaknesses. New, complementary or alternative tools that have been developed in the project feeding this paper are presented, which can contribute to the scientific community and stakeholders of the topic. The main result of our research is a structured and practical guide (a ToolBox Manual) that can support practitioners and researchers wishing to characterize and classify large rivers, based on the River Styles Framework. The main contribution is that this set of ideas, solutions, and tools, makes this type of exercise significantly more transparent and at the same time much less subjective. Moreover, the procedure is applicable to large systems and does not require more information than that generally available also in developing or emerging countries.

1. Introduction

Today, remote sensing offers unprecedented opportunities to analyze rivers, particularly in large basins. Piégay et al. [1] present a very complete, articulated, and up-to-date panorama of the techniques, recent experiences, and emerging research challenges that use or are related to remotely sensed data for geomorphic river characterization and management and, particularly, to investigate past, present, and future fluvial corridor conditions and processes. Most of these very refined and complex techniques and experiences relate to the identification of geomorphic units (e.g., main channel, bars, islands, etc.) or the measure of some features (e.g., topographic, vegetation structure); others are focused on the search for (cause–effect) relationships between natural or anthropogenic factors and the geomorphic response of the river, including its time evolution (e.g., widening, narrowing; incision/aggradation); still others are dedicated to monitoring river behavior, either through periodic, systematically repeated surveys or real time observation. It can be argued that there is no better way to describe how a river actually works than observing its actual dynamic behavior along a significant time span and across different spatial scales. However, a basic, fundamental understanding can be obtained by looking at the main, present characteristics of the river at different scales, simultaneously, and recognizing how different assemblages of geomorphic units occur in relation to forcing factors (“controls”) such as geology, topography, climate, and human interventions. This essentially static analysis, when integrated with information concerning the historical evolution of a river, provides a very sound basis for understanding its character and behavior from a geomorphological perspective. This is the essence of the River Styles Framework (RS onwards) developed by Brierley and Fryirs [2]. Of course, information generated with this framework can be corroborated, refined, and integrated with further, more advanced and complex analyses, where tools like those presented by Piégay et al. [1] can be extremely valuable.

The River Styles approach, however, requires a quite demanding analysis effort, particularly when dealing with large basins. In addition, it incorporates a high dose of subjectivity as it involves several expert-based judgment inputs. This latter point is a delicate one, as computer-based tools and remotely sensed data evolve and somehow challenge the role of the human expert (geomorphologist), which is hence being continuously re-discussed [3,4].

The aim of this paper is widening the audience of river practitioners, managers, and researchers that would like to apply the RS Framework as potential beneficiaries, but which may be reluctant to do because of the practical difficulties that are typically encountered, particularly when trying to apply it to large river basins. We contribute to this aim, on the one hand, by defining a more structured, objective, and straightforward procedure, in which the role of the geomorphology experts is more clearly identified and circumscribed, and, on the other hand, by suggesting concrete computer-based tools options to support the required tasks and, even, automate a large part of the exercise. While Fryirs et al. [4] offer a broad spectrum view of available tools and their possible role in the RS analysis (and more), which is certainly stimulating, we concentrate on the practitioners’ point of view, as we assume they are keen to find as simple and concrete guidelines as possible, to actually perform the multiple tasks.

It is important to point out that the RS Framework involves a number of stages and steps, of which we only consider here “Stage I, steps 1 and 2”, i.e., the “characterization-classification” exercise; in what follows, we refer to these as the “RS analysis”. The findings presented in this paper are products of a project conducted by some of us for The Nature Conservancy, in Colombia; we are referring to it in the following as the “GeoMagda project”. It is worth noting that these findings have been articulated in a number of more or less simultaneous articles and, as such, unavoidably we need to refer to them.

Several important research groups are contributing to the plethora of methodologies and tools aiming at supporting geomorphic river analysis. With no pretention of exhaustiveness, we identify the following main ones:

◦

The Riverscapes consortium https://www.riverscapes.xyz/, a collaboration between fluvial researchers in the United States of America, United Kingdom, New Zealand, Australia, Italy, France, and Spain, particularly with its ToolBox RCAT (http://rcat.riverscapes.xyz/), which includes: the Valley Bottom Extraction Tool (V-BET), the Riparian Vegetation Departure (RVD) tool, the Riparian Condition Assessment (RCA) tool, and the Riparian Recovery Potential (RRP) tool.

◦

The Valley Bottom Confinement Tool (VBCT), developed and based on studies by Fryirs et al. [4,5] and O’Brien et al. [6] (http://confinement.riverscapes.xyz/).

◦

The TerEx ToolBox (terrace and floodplain mapping), based on the work of Clubb et al. [7] and Stout and Belmont [8] (https://qcnr.usu.edu/labs/belmont_lab/resources).

◦

The Fluvial Corridor ToolBox (FCT), by the group led by Hervé Piégay at the University of Lyon (France); currently, they are developing a version within the QGIS 3.x interface compiling scripts written in Python (https://github.com/tramebleue/fct).

◦

The River Meandering Migration Software (RVRMeander; http://www.rvrmeander.org/) of the University of Illinois, University of Pittsburgh, and USDA(U.S. Department Of Agriculture), with applications in particular to the Amazon River by Jorge Abad.

◦

The Geomorphic Unit Tool (GUT), which delineates instream geomorphic units (GUs) from topography using a three-tiered hierarchical classification adapted from Wheaton et al. [9] (http://gut.riverscapes.xyz/).

To this we can add the REFORM (REstoring rivers FOR effective catchment Management) project [10] presenting a host of methods concerning the interpretation of remote sensing data (e.g., [3,11,12,13,14,15]), as well as the already mentioned Fryirs et al. [4] and Piégay et al. [1].

The use of most of these tools, however, is typically quite complicated and delicate. They generally require a Python installation that often generates difficulties. They always require a human refinement of their outputs, which introduces a hidden degree of arbitrariness. Some (such as RCAT) are explicitly designed to be fed with the standard information available in the USA. Others (e.g. the RVRMeander) are explicitly aimed at meandering rivers. Several are in principle interchangeable because they address the same task, although with some differences. No single one, however, embeds all of the functionalities needed to carry out a River Styles analysis. Finally, most of them are simply not applicable in cases characterized by scarce information, including low resolution geographic data (particularly Digital Elevation Models-DEM).

The most applicable tool to support RS analysis, amongst those listed, is certainly the Fluvial Corridor ToolBox. However, its actual usefulness for performing RS analysis is definitely much reduced, for several reasons, with respect to the potentialities that were originally been announced in the FCT Project. On the one hand, the sequence of sub-tools that appears when installing it (within ArcGIS) does not correspond to the actual sequence of use, a quite confusing fact. Additionally, several steps and functionalities are not relevant for RS analysis and thus must be skipped, while, conversely, other key functionalities are not supported (e.g., confinement or RS classification). Another issue is that the computation of sinuosity within the FCT is based on the polyline connecting inflection points of the active bed centerline; in some cases, this may result in a polyline external to the valley bottom (VB), which should not happen. Finally, the key Hubert test for identifying homogenous reaches does not work with non-cardinal indicators, which are nonetheless often involved in RS analysis. Moreover, several “tricks” need to be known; to a novel user, many of these may appear to be unsolvable snags, leading him/her to quit.

This paper addresses mainly the community of practitioners by providing a guide through the different elaborations needed when conducting a RS exercise. Its main contribution in this sense is providing a procedure that allows to link all the required steps in a structured, operational form. In such a way, on the one side, the role for expert judgment becomes better defined, so that the outputs are structurally more objective; on the other side, the procedure makes RS characterization–classification viable also for large rivers, by taking advantage of computer-aided tools. The paper may also be of interest to the scientific community, as it addresses some challenging conceptual problems concerning the holistic synthesis of reductionist information that lie at the core of the automated analysis of complex information. We are aware, anyway, that the main added value of this paper can be appreciated by those who know and have adopted (or would like to adopt) the River Styles Framework.

In what follows, we first present the essence of the proposed procedure (details are provided in Appendix A); then, some of the newly developed tools are presented, together with demonstrative applications to the Magdalena River case study (Colombia); finally, conclusions are drawn.

2. Systematizing a River Styles Inspired Characterization–Classification Exercise

Figure 1 presents a summary view of the steps we consider necessary to conduct a RS analysis (Stage I, steps 1 and 2 of the RS Framework), while

Table A1. Steps operationally needed to carry out a full River Styles (Stage I, Steps 1 and 2) analysis and corresponding tools/procedures proposed for doing so. Bold italic type items are those we developed “ex-novo” and which are fully described in the ToolBox Manual; two of these—“Confinement tool” and “Reductionist → Holistic algorithm”—have been described in the text above (Section 3.1). Notice that in our application to the Magdalena River Basin, described in [16,17], we adopted a set of official digital ArcGIS shapefiles that allowed us to skip some of the steps related to the delineation of geomorphic elements and units (but not all), but at the same time raised some issues mainly due to resolution and lack of synchronicity with the available Digital Elevation Model (DEM) from the Shuttle Radar Topography Mission (SRTM) 30 m filled.

	Object	Algorithm	Tool/Procedure	Observations
Methodology	Procedural Tree and definition of attributes	Choose structure suited for the case at hand conserving a river basin (or even better, a regional) perspective	Expert judgement	This is where expert judgment is mainly involved
Instrumental	DEM clipped to the area of interest	Select and extract zone of interest	Manual: ArcGIS TOOL CLIP (Data Management) according to a rectangle mask or EXTRACT by MASK, given a polygon (i.e., the river basin)
Instrumental	Landforms (terraces, fans…)	Expert judgement (from aerial photographs, field surveys, …)	Manual	With a high-resolution DEM some degree of automation is possible (e.g., TerEX by Stout and Belmont, [8])
Instrumental	GIS reference points	Locate points along the axis of VB (see below) at constant distance	Standard GIS procedures (creating points along a line)	These are very useful to orient the analyst while quickly switching from GIS to Google Earth
Instrumental	Geographical reference points	Expert judgement (towns, known sites, appropriate density, closeness to river, manageable number, etc.)	Manual construction of a point GIS Shapefile (SHP)	Sites useful to guide descriptions and to show in schematics and maps
Instrumental	Base segmentation points	Location where there is an evident change of VB width, catchment area (tributary); break in slope; discontinuities (dams)	Manual standard GIS procedures. Possibly use FCT: Disaggregation; Metrics (width, etc.), Statistics (Hubert’s test) and Metrics (Discontinuities)	Can be used as a preliminary indication to define reaches (e.g., VB width well identifies macro reaches)
Instrumental	Type of VB margin	Identification of infrastructure and landform mapping	Manual by standard GIS procedures	First the effective valley bottom must be defined (see below)
Instrumental	Contact points river-valley margin	Identify key points where the envelope of active riverbed (see below) touches the margin of the effective VB. Geological, planform, and anthropogenic contacts should be considered	Manual, expert judgment procedure on GIS supported by DEM info + standard construction of a GIS point SHP (semi-automatic GIS algorithms are an alternative)	Useful to assess degree of "constraint" and sediment sources, particularly in laterally unconfined (but constrained) reaches. Active channels’ envelope needs to be determined first (see below)
Instrumental	Contact segments between active riverbed and VB margins	Create a thin buffer around the envelope of the active riverbed; then intersect with the effective VB (see below)	Manual standard ArcGIS procedures (buffer, intersection)	It is a basic information for the confinement attribute. Active channels’ envelope must be determined first (see below)
Instrumental	Sediment sources	Identify where there is: a contact point with landslides, fans or terraces; or meanders with external VB elevation higher than internal (so that the erosion-sedimentation balance is positive); or a tributary input	Manual (expert based) or semi-automated GIS procedures	The expert judgment can be synthesized in a logical tree that can then be implemented (even in an Excel spreadsheet) to carry out this task
Instrumental	Constrictions	Identify where:-VB is narrow: (V-W local)/W < α-the reach is single thread and not excessively "spread"; i.e.,: [ W–w]/w < β -The envelope significantly reduces with respect to the upstream average: [W–W ]/W > γ. Where: CA: active river bed (set of active channels); V: VB width (and V: moving average); W: width of envelope (W: moving average); w: width of active river bed; and: α,β,γ: parameters.	Excel spreadsheet (GeoMagdaToolBoxManual: "Esqueleto_magda")	A constriction is not a contact point, but a short reach where the VB narrows down significantly, at the local scale, as compared with upstream and downstream. Constrictions are important because they can induce changes in hydraulic and sediment behavior and, consequently, in morphology
Instrumental	Segmented VB	Spatial discretization of the VB according to regular reaches along its centerline with a relatively fine step (e.g., 500 m for the Magdalena) and sequenced (i.e., spatially ordered)	FCT: Disaggregation; Sequencing	This is a very important step that supports several subsequent steps. The key output is the “Rank_DGO” field of the obtained table.
Geomorphic elements	Active channel “just water”	Satellite images interpretation. For Magdalena: use of official SHP IGAC 100 mil “Drenaje doble” already available; separation of branches not completely connected and tributaries	-ArcGIS procedure: calculation of indices like MNDWI (the Modified Normalized Difference Water Index) [30,31], selection of pixels over threshold, transformation to polygon, cleaning-Manual, expert judgment in ArcGIS	It is important to use images referring to high waters, but actual overbank flooding should be avoided. Automatic algorithms suffer from presence of clouds
Geomorphic elements	Active riverbed (CA) (collection of active channels completely connected, including bars, but not islands)	Add mid channel and bank attached bars to the "active channel just water” polygon. NOTICE: The procedure presented here refers to the specific information available for the Magdalena, where an official SHP IGAC 1:100,000 was available for bars ("bancos de arena") but did not discriminate amongst main river and tributaries. It must be noted that this SHP presents some inconsistencies; it was probably built based on different, asynchronous images, possibly by different operators.	ArcGIS procedure: -draft (still without bars) envelope of active channels (see next step); -buffer of this envelope of a max width of lateral bar; -CLIP the original bars SHP (IGAC 2017 100 mil «bancos de arena») with this envelope; -manual check to eliminate bars not belonging to the river of interest; -MERGE with the active river bed; DISSOLVE	Care must be put in not including bars belonging to tributaries; there is a degree of ambiguity. A more direct procedure consists of identifying bars directly from satellite images as bare surfaces, based for instance on the SWIR (Short-wavelength infrared) band (band 7 in OLI (The Operational Land Imager), TIRS (Thermal Infrared) and TM (Thematic Mapper) Landsat missions), as shown in Monegaglia et al. [32]
Geomorphic elements	Envelope of active riverbed (collection of connected active channels)	Transform the polygon into a polyline SHP and eliminate all lines different from the external boundary	ArcGIS procedure:-Tool «Feature to Line»; -manually eliminate internal small polygons; -retransform into a polygon with «Feature to Polygon»; -finally, possibly use «Smooth Polygon» (Cartography)	This is not a "convex hull", but the polygon ideally obtained by merging the polygons of all active channels fully connected
Geomorphic elements	Islands	Eliminate the active riverbed (including bars) from the envelope of active riverbed	ArcGIS: ERASE	Alternatively, interpretation of satellite images with appropriate indices (e.g., NDWI, Carlson and Ripley, [33]) can be carried out. However, the procedure suggested here ensures consistency
Geomorphic elements	Bank-attached bars	Eliminate the envelope of active riverbed "just water" (which, by construction, includes islands and mid channel bars) from the envelope of active riverbed (CA)	ArcGIS: ERASE	Alternative methods exist; for instance, the GUT tool within the River Scape Community ([34,35,36]), but they just work with high resolution imagery and altimetry. UAV surveys can provide viable alternative at moderate scales of analysis (e.g., Casado et al. [12] or for a review: Carrivick and Smith, [3])
Geomorphic elements	Mid-channel bars	Eliminate the active riverbed "just water" and island sand bank attached bars from the envelope of active riverbed	ArcGIS: ERASE triple	Same observation as above
Geomorphic elements	Axis of active riverbed	Centerline of the envelope of active riverbed	-FCT (Tool 1.3 Centerline), or ToolBox Manual (Módulo Básico 2): tool based on ArcGIS/ArcScan which provides a more smoothed output-PyRIS proposed by Monegaglia et al. [32]: the tool identifies the axis of the main channel even in anabranching rivers, defined either as the widest or the longest	For RS application, a centerline is sufficient and perhaps preferable because in anabranching rivers, selecting a particular branch may be less representative.
Geomorphic elements	Natural valley bottom (VB)	Identify the portion of flat land adjacent to the river and not higher than a given threshold. Build to this aim various DEM of Differences (DoD) with various thresholds and check consistency with land-use, envelope of contemporary flooded areas (e.g., from the Dartmouth Flood Observatory database: https://www.dartmouth.edu/~floods/, specifically the derived product “Annual flooded areas”), geomorphic map, field evidences	FCT: Spatial Components; 1.2 valley bottom. See however various details explained in the ToolBox Manual, "Modulo VB". Other tools exist like for instance VBET [37], but with no particular advantages	Even with a perfect DEM, a manual correction is needed because of several reasons, e.g., the DEM may have captured the elevation of water, not of riverbed and hence the DoD is at least ambiguous.
Geomorphic elements	Effective VB	Eliminate from Natural VB portions separated by works that may interrupt channel migration or flooding (e.g., roads, railways, rip-rap protections)	Manual procedure in ArcGIS. In the ToolBox Manual, "Modulo Básico 8", a systematic procedure is proposed	-
Geomorphic elements	Axis of effective VB	Centerline of the envelope of effective valley bottom	Totally analogous to "Axis of active riverbed"	-
RS Attributes	Sediments (bed material texture)	Construction of a point SHP in each sampling or survey station and transformation to a line SHP (based on the Axis of Active Channel) assuming some interpolation criterion	-Physically sample or adopt some imagery interpretation tool (e.g., [38,39]) -Interpolation: assume that the sampled texture holds until next station, possibly with expert-judgment adjustments; or: use an algorithm to translate reductionist information into a holistic synthesis like the "HOLISTICO binario" Excel Tool	-
RS Attributes	Planform	Identify reaches with a recognizable combined pattern of number of channels, distance, width, sinuosity and number, size and position of islands and bars.	-Expert-based identification based on satellite imagery (e.g., Google Earth) and available polygons of active riverbed. Supported by ToolBox Manual: “Modulo FCT 3.4 Morphometry” for sinuosity and “Modulo sinuosidad nueva”; “Modulo MonoMulti canales” for single/multi thread (see this paper). -Alternatively, use the automated Excel tool presented in Nardini et al. [16]	a-Single thread: rectilineous, sinuous, meandering, irregular, tortuous. b-Transitional: wandering; alternate bars; “swallowing” (i.e., still single thread and mainly straight or sinuous, with local significant widening, including one or more islands. c- Multi thread: braided; island braided; anabranching avulsive; anastomosed; combinations.
RS Attributes	Streamflow character	Identify reaches with a recognizable pattern of current threads, turbulence, areas of still water. (for rivers where the riverbed is visible at low flow, a similar attribute would be “bed morphology”)	-Expert judgement or - Other Tools like GUT (see above) -interpolate with reductionist-holistic “HOLISTIC categorical” Excel Tool	Google Earth imagery brings in ambiguities because of varying quality, date and time of the images.
RS Attributes	Presence of natural levees	Identification of elements from aerial or satellite imagery (e.g., Google Earth) and DEM; holistic synthesis of the reductionist information into meaningful stretches; translation reporting into a polyline along active riverbed axis with binary attribute (presence/absence)	Expert judgement + info from existing reports → construction of a polygon SHP Then: GIS intersection of such SHP with the segmented effective VB (from FCT); translation into a binary attribute (presence/absence) along the discretized VB axis (GIS procedure); finally, translation into a polyline (presence/absence) along the VB axis with the “ HOLISTICO binario” Excel Tool (see above); possibly, transposition into a polyline along the axis of the active river bed with the Excel Tool “De_VB_a_CA”	For the first step of identification, alternatively, GUT Tool or similar can possibly be used. The transposition from the discretized VB axis into the axis of the active river bed is not trivial because, owing to their different location in the floodplain, a different number of segments from discretization may occur. All this is explained in detail in the ToolBox Manual (Modulo Basico 5 and 6)
RS A.	Presence of wetlands	analogous	analogous	analogous
RS A.	Presence of paleo channels	analogous	analogous	analogous
RS A.	Presence of oxbows	analogous	analogous	analogous
RS A.	Presence of ridge & swale topography	analogous	analogous	analogous
RS A.	Presence of bank-attached bars	analogous	analogous	analogous
RS A.	Presence of mid-channel bars	analogous	analogous	analogous
RS A.	Presence of islands	analogous	analogous	analogous
RS Attributes	Confinement	Identification of contact lengths (see above) and computation of a statistic indicator over the reaches of the selected segmentation	Manual ArcGIS procedure as detailed in the ToolBox Manual Modulo 3 bis (this paper)	Confinement is a statistical indicator defined over an ordinal scale, assessed for each reach, with values: “confined”, “partly confined”, “laterally unconfined”
RS Attributes	Cause of confinement	For each confined or partly confined reach, identifies which is the prevailing cause amongst: “valley”, “planform”, “works” (anthropogenic)	Manual ArcGIS procedure as detailed in the ToolBox Manual Modulo 3 bis (this paper)	In each reach, several causes may intervene; the prevailing one in terms of lengths involved (adding right and left bank) is considered
RS A.	Constraints	Expert judgement based on the contact points already identified	Manual on ArcGIS
Synthesis	Reaches	According to Nardini et al. [16], we identify independent segments according to non-extensive attributes: Planform; Bed material texture; Streamflow character, finding the “least common, spatial denominator”	(see Nardini et al. [18])	The reaches are the basis for the Geo database and to support visualization of attributes
Synthesis	GeoDB (polyline GIS SHP segmented according to reaches)	Assign to each reach a value for each one of the relevant attributes. For those directly involved in the definition of reaches, just select the appropriate value associated with the original segmentation; for the remaining ones (belonging to the Procedural Tree), compute a statistic, ordinal indicator (e.g., for binary cases like local presence/absence of a geomorphic unit: 1: “absent”, 2: “occasional”, 3: “frequent”, 4: “prevailing”, 5: "complete”)	For: Cause of confinement: see Excel Tool “CONFINA” in the ToolBox Manual. For binary attributes it is the % of reach length where the attributes presents a “YES” value: see Excel Tool ESTADIS_Tramos	According to our Procedural Tree, the attributes requiring a statistical computation are:-confinement-cause of confinement-geomorphic units within the VB (levees, oxbow, ridge and swale, paleoCH, wetlands)-geomorphic units within the channels (bank attached and mid channel bars, islands; streamflow character)
Synthesis	Skeleton (finely discretized, polyline GIS SHP)	Discretize the VB axis in short segments (e.g., 500 m for the 1400 km long VB axis of the Magdalena River); compute a metric for each attribute relevant to generate a meaningful synoptic scheme (distances, elevation, contributing catchment area; active river bed width, VB width, valley slope); build an area- bankfull flowrate relationship (Qb = Qb(A)) and determine total and specific streampower	FCT: 3. Metrics (in each segment, the Tool computes several values, e.g., every 50 m within a 500 m segment, and offers sub-statistics)	Skeleton is better constructed along the axis of the valley bottom in order to represent planform sinuosity. In the Magdalena exercise, owing to the low DEM resolution, FCT did not work to determine the contributing catchment area; a manual ArcGIS procedure was adopted instead and then data were processed within the Excel “Tool ESQUELETO”
Synthesis	Synoptic scheme (profile and plan view)	Build a sketch in Excel (both profile and plan view) by using the Skeleton information and reference points SHP.	A suitable filter (e.g., moving average plus not decreasing condition) must always be applied to the original elevation data from DEM, particularly when referring to the VB axis because this is an idealized line that often crosses planforms or hills (an example is offered in the Excel Tool ESQUELETO)	There must be one for each river considered in the basin (in our preliminary exercise Magdalena and Cauca rivers). It is appropriate to include in the main river sketch (Magdalena) the long profile for all tributaries
Synthesis	Macro scale map	Map the relevant information (possibly available from existing cartography): landscape units, faults involving the rivers; active channels; reference points; effective VB; macro scale confinement; planform (or at least single-multi thread segments); cause of confinement (VB margins); constrictions	Standard ArcGIS mapping procedures	All relevant rivers are presented in a same map. This is split into a minimum number of whole basin maps because merging too much information is confusing; so, the clearest, but most synthetic combination must be sought
Classification	River Styles	Identify the reaches with the same combination of values of main attributes, according to the Procedural Tree adopted (e.g., planform, bed material texture, confinement, cause of confinement); do the same with all attributes.	Python/ArcGIS Tool (ToolBox Manual): see ToolBox Manual and Nardini et al. [16,17]	Owing to ArcGIS limitations, currently it supports only 7 attributes at a time; therefore, the process might need to be split in more than just two cycles. Actually, we needed three for the Magdalena:I cycle: main ones → Label Main II cycle: Label Main + 5 geomorphic units within the VB → Label II. III cycle: Label_II + 4 geomorphic units within the channel

in Appendix A presents a much more detailed picture. For each step, we indicate the tool/procedure suggested; in some cases, several options are proposed, depending on the original information available. The list does not claim to be exhaustive (see also [1,3,4,11,12,13,14,15] for additional options); rather, its purpose is to give an idea of the complexity of a full RS analysis exercise, and to locate within the process the additional tools we explicitly developed, which are described in detail in the ToolBox Manual [16].

Expert judgment is involved on one hand in defining criteria to be adopted within the various steps (e.g., defining the Procedural Tree; or defining how to distinguish a planform type from another) and, on the other, after the procedure has been applied, in interpreting river behavior. This exercise must be accompanied by an analysis of the river’s historical behavior and recent adjustments (not included herein).

As anticipated, the added value of this procedure can be fully appreciated by those who know and have used (or would like to adopt) the River Styles Framework. The innovations are evident only when analyzing the details, but consist mainly of:

(a)

Modifying the Procedural Tree, by “locating first what weighs more”, in such a way that a hierarchical classification can be obtained that can manage all sub-types generated in particular by the assemblage of the different geomorphological units (which can give rise to a large number of sub styles).

(b)

Identifying objectively the (variable length) reaches according to the spatial intersection of only the set of intensive-type attributes of the Procedural Tree (i.e., those not depending on the extension of the length of the considered segment) and assessing all the extensive-type attributes as spatial statistics over such reaches. This idea is further explained in [16,17,18].

(c)

Relying on a spatial database (called “skeleton”) obtained by segmenting the valley bottom in sub-segments along its axis, and determining a value for each relevant attribute in each sub-segment.

(d)

Specifying a set of tools that ensure the ability to also deal with large river systems. The proposed tools (Appendix A) are a particular choice out of a virtually infinite (or at least quite large) number of alternatives; but having defined at least a tool set that does cover all the demanding steps (elaborations of extensive information and computational steps) brings in an added value as, without them, several snags might prevent a candidate user from applying the RS approach.

(e)

Proposing a sequence of operations that covers the whole exercise, ensuring a linkage amongst steps and avoiding redundancies.

Concerning aspect (d), the ToolBox Manual developed within the GeoMagda project provides quite detailed instructions on how to solve a number of difficulties; this might be superfluous to the trained GIS expert, but precious to a broader arena of users. It is worth stressing that the tools we propose require just a basic expertise in computer use, particularly GIS and Excel; this has a negative facet, as advantage is not taken of more powerful possibilities, but it also has a positive facet, as broadens the arena of candidate users and to improvements.

It is worth noting that the first main product of this paper is Table A1 in Appendix A, which presents the full procedure we propose.

In the following section, we address just three “tools” amongst those developed and included in Table A1, in order to give an idea of the type of contribution this paper (and the underlying GeoMagda project) brings in. To illustrate some of the concepts and type of results we utilize the Magdalena case study developed in [16,17]. What follows should not be interpreted, however, as the application to the Magdalena of the full procedure; rather, the focus here is using the Magdalena data just to illustrate how these “tools” work (any other data would serve this purpose). A synthesis of the Magdalena case study is available in Nardini et al. [17], where we show the essential outputs obtainable by applying the procedure described herein; details can be found in [16,17]. We warmly suggest that this paper be read together with the case study (and also with Nardini et al. [18] when it be published) because they are indeed highly complementary. From a methodological point of view, Nardini et al. [17] present the essential contents of a River Styles characterization–classification exercise, while here we present the organized collection of technical steps involved.

It is important that the reader not get lost in the following technicalities: the essential message we would like to convey is that the procedure we have developed addresses several points similar to those discussed below, proposing conceptual and operational solutions to overcome any issues, while simultaneously ensuring that all steps are linked in a harmonic fashion.

3. New Tools/Procedures

3.1. Confinement

Confinement, as adopted in the RS Framework, expresses the extent to which the river abuts either valley bottom (VB) margin. It is a very meaningful attribute because this condition strongly characterizes the dynamic behavior of the river itself. It should not be confused with “entrenchment”, i.e., the ratio between active channel width and VB width [19].

To assess confinement, the objective, relevant base information is the set of contacts between the active channel (or, more generally, the envelope of active channels) and the VB margins. However, to meaningfully characterize a reach in terms of confinement, one has to rely on a statistical, categorical indicator (% of either bank’s contact lengths over the corresponding axis length) according to some given convention (according to Brierley and Fryirs, [2]: laterally “unconfined” when total contact is less than 10%; “confined”, when is over 85%; “partly confined”, when in between) over a given reach (product of a segmentation). Indeed, the result strongly depends on the segmentation adopted, as shown in Figure 2.

The main novelty of the algorithm we propose to define confinement lies in determining it after reaches are defined. In fact, this position overcomes an intrinsic logical problem of the traditional approach: in the RS framework, reaches displaying a recognizable and meaningful pattern are identified as the output of a process based on assessing the set of attributes included within the RS Procedural Tree; but the first of such attributes is indeed confinement, which requires—as we have shown—that segments (i.e., reaches) be pre-defined.

Additionally, our algorithm determines an explicit statistical confinement indicator as well as another important attribute which is the cause of confinement, once the polylines of the VB margins separated in categories (i.e., valley margin, planform, and infrastructure) are provided. Furthermore, it brings a clear advantage from a practical point of view: it is accessible to anyone able to use ArcGIS (or QGIS) and Excel.

Finally, as confinement is determined for any given segmentation, it is possible and meaningful to provide the output at different scales; particularly:

(i)

reaches (as usually done in the RS);

(ii)

basin scale (macro): indeed, it is possible to adopt a very simple and effective criterion to identify macro segments, e.g., where a significant change of VB width (V) occurs (a fact that can be detected by the Hubert Test in the FCT or simply...visually) and the resulting macro confinement is very meaningful (Figure 3).

As O’Brien et al. [6] have already presented an algorithm, with a corresponding GIS Python software tool, to determine the confinement indicator within the RS context, it is important to recall it and point out the differences with our approach. That algorithm goes as follows:

◦

it determines the contacts—both margins—between active channel and “inflated” (buffer) VB margin through a GIS intersection;

◦

these contacts are projected onto the reference axis with a binary value (0 for no contact; 1 for contact at either bank);

◦

the axis is discretized in “short-segments” (identified by k) of a pre-fixed length;

◦

the algorithm goes on by computing locally (for each slice) the statistical indicator c(k) as the total length of the projected contacts (either banks) falling in it, over the segment length;

◦

finally, it assumes that when the indicator c(k) is 0.50 < c(k) < 0.85 the cause of confinement is the valley margin, while when it is 0.10 < c(k) < 0.50 the cause is a planform geomorphic unit.

O’Brien’s algorithm implements a peculiar way to obtain a holistic synthesis, starting from reductionist analytical information associated with short river segments. We found that this may lead to inconsistencies (as shown in Figure 4) or even mistakes (possibly due to our inappropriate use?), when the river is anastomosed with a very wide envelope of the active channels (“AC envelope” in what follows) (see Figure 5). Moreover, the assumption concerning the cause of confinement is evidently arbitrary and, furthermore, cannot include the frequent case of confinement by infrastructure (e.g., railways, roads, levees). Finally, it does not compute an explicit statistical indicator: the classification is obtained only by assigning a color scale.

In general, the algorithm produces a somehow hysterical sequence of different confinement categories (in this case: mostly laterally unconfined with two shorter segments partly confined) and is not able to capture a more holistic character (indeed, we would classify the whole stretch as partly confined). Notice in Figure 4 that while the output of the algorithm is reported along the reference axis adopted (colored polyline axis: dark green: unconfined; light green: partly confined planform controlled; orange: partly confined valley margin controlled; red: confined (not detected in this figure along the axis)), to assess confinement one has to look at the blue polygon representing the envelope of active channels (which is wide or very wide where the river is anabranching or anastomosed), and identify sectors where it touches either margin of the VB polygon (grey area). Only the red contacts count because even though other lines seem in contact, visually, from a geographical point of view (GIS) they are not: sometimes, indeed, the difference amounts to a few meters only, but it is still larger than the tolerance adopted in the buffer used by the algorithm.

The algorithm we developed in GeoMagda overcomes such difficulties. As already stated, it assumes that the reaches over which the confinement statistical indicator is to be computed are already defined. This is the basic difference with O’Brien’s algorithm.

Other differences, which constitute additional advantages, are as follows:

◦

it is accessible to any user of ArcGIS (or QGIS) and Excel;

◦

it computes an explicit statistical confinement indicator;

◦

it determines the cause of confinement, if provided with the polyline of the VB margins separated in categories: valley margin, planform, or infrastructure.

Our algorithm undoubtedly suffers from some limitations: It is designed for isolated rivers (not a river network); the current version is based on the number of discretized short-segments with a contact (the “numerosity”), not on the actual contacts length; and it is not automated, which means that it requires a rather large number of manual (but not expert based) steps that, as a whole, can be quite cumbersome. However, with a sufficiently fine discretization, there is no actual gain in considering the length of the short segments rather than their numerosity; and automatizing such a procedure is perfectly possible.

In Appendix B we provide details of the algorithm, so that it can be replicated step by step; Figure 6 provides a general idea of the first part, related only to confinement.

3.2. A Logical Heuristic “Reductionist→ Holistic Algorithm” for Mono-Dimensional Variables

We address here the need to obtain a holistic synthesis from reductionist information. The expression “Reductionist → Holistic” refers to the attempt to translate information related to an attribute at the local level (e.g., the width of the active channel or the presence/absence of islands within a short segment of the river) into information that synthesizes, in general, multiple attributes and with a view that spans well beyond the single short-segment. In the particular mono-dimensional applications presented below, the reductionist information is the local presence/absence value (binary case) or else the specific value within a finite set of possibilities (categorical case), while the holistic synthesis is a judgment obtained by looking at a sequence of such reductionist values along the river, identifying homogenous stretches. To this aim, we developed a logical heuristic algorithm, both for the binary and thee categorical cases, whose results are nicely consistent with expert judgment output. We adopt the term “logical heuristic” to stress the fact that the algorithm is based on explicit reasoning (logics), but with some steps that not necessarily lead to an optimum, although still achieving the purpose (at least in our application).

3.2.1. Reductionist → Holistic Algorithms: The Literature Status

For the mono-dimensional problem (only one attribute), a statistical approach seems to be the best option. Parker et al. [20] provide an interesting literature synthesis of statistical algorithms borrowed from geological analysis [21], recalling Davis’ classification scheme [22] in (i) “local boundary hunting” and (ii) “global zonation”. The former searches for the boundary between two distinct levels of the considered character; e.g., Webster’s algorithm [23] and its variants identify the location of “breaks” by sliding a moving window until the difference of variance before and after its midpoint is maximized, while variance is kept at a minimum within each one of the two halves. it’s the performance depends on the width of the moving window, a variable which is selected somewhat arbitrarily. The second approach, “global zonation”, looks at the whole data set at a glance, trying to identify those subsets of points that are internally most uniform. For instance, Bohling et al.’s [24] cluster algorithm identifies first that location characterized by the minimum change of any possible couple of adjacent points within the dataset. These two adjacent zones with the smallest difference are then combined into one combined zone, which is treated as a single object, characterized by the mean value of the points within the zone, and considered in the next iteration. The process continues, with more and more zones being joined based on their similarity to each other. This method begins with zero within-zone variance, but each iteration results in an increase, until the proportion of total variance explained by the zonation falls below a specified threshold. Other “global zonation” methods, however, like Gill, [25] and Hawkins and Merriam, [26], are better suited for the identification of river reach boundaries because they statistically minimize within reach variation and maximize between reach differences, and are less influenced by local inconsistencies in the data sequence. Gill’s method explicitly seeks to identify the zonation that minimizes the variance within each zone (reach) and maximizes the difference between the zones (reaches). Hubert [27] introduced an algorithm which works in a similar manner; although it is possibly the best performing one for this task, its outputs are not exempt from questioning (Figure 7). Similar conclusions have been found by other researchers (e.g., Martinez-Fernandez et al. [28]).

Furthermore, it must be recognized that all of these approaches have been designed thinking of quantitative, continuous (scalar) variables; when binary variables are concerned instead, they conceptually and practically fail. This is particularly true for Hubert’s test, which is considered the key tool to define separate reaches within the Fluvial Corridor ToolBox [29]. This is one of the reasons why we explored a different approach by developing logical heuristic algorithms.

For evaluating the performance of algorithms, Parker et al. [20] assume that the higher the proportion of variability R is explained by the reach boundaries, for a given number of reaches, the stronger the performance, and hence the zonation algorithm is best suited to defining functional river reaches. Operationally (what follows has been slightly reformulated and corrected from Parker et al. [20]) R is equal to Equation (1):

R = σ²_zone/(σ²_zone + σ²_stretch)

(1)

where σ²_zone is a measure of the variance of the reach means about the grand mean μ of the whole sequence, defined as follows Equation (2):

σ²_zone = [∑_j=1…m (μ_j − μ)²]/(m − 1)

(2)

where σ²_stretch, j, m and μ_j are defined as follows Equation (3):

σ²_stretch = ∑_j=1…m[ ∑_i=1…nj (x_ij − μ_j)²]/[N − m]

(3)

where ${\mathrm{x}}_{\mathrm{ij}}$ is the $\mathrm{i}$th point within zone $\mathrm{j}$, μ_j is the mean of the $\mathrm{j}$th zone, ${\mathrm{n}}_{\mathrm{j}}$ is the number of points in the $\mathrm{j}$th zone, N = (∑j=1…n_j) is the total number of points considered amongst all zones, and $\mathrm{m}$ is the number of zones.

The closer the performance index R gets to 1, the better the output is, because this means that the intra-cluster variances countless (i.e., they are smaller because cluster data are more similar to cluster mean), while the difference amongst clusters (expressed as the variance of the clusters means) increases and becomes the dominant component of the differentiation. Thus, one has clusters internally more homogeneous and more diverse amongst them.

Although very reasonable, this criterion is not suitable for our case. Indeed, as is apparent from what follows, our algorithm tries to “give life” to all cases of “presence” as soon as they form a sequence sufficiently long (actually longer than L_M), without leaving “holes” longer than L_M, being this a significant length parameter to be set. So, the performance must be assessed on this basis. At the moment, we just considered ourselves satisfied with the visual appreciation.

3.2.2. Reductionist→ Holistic Algorithm: The Binary Case

In this case, the variable of interest is assessed at a number of points along the river in terms of presence/absence (usually corresponding to a relatively dense discretization as presented in Figure A1–Appendix A). Hence, the algorithm presented here operates on a binary variable (0,1) that maps the original variable over a fine discretization of the river axis into segments. We assume that the character of the river is assessed over a stretch of meaningful length, in no case shorter than a stretch with a “minimum significant length” L_S. The philosophy interpreted by the proposed algorithm is hence in a few words: “identify those stretches, not shorter than L_S, to which a predominant value can be assigned”.

The algorithm operates as follows:

◦

choose a significant minimum river length LS, basically according to the usual criterion of LS = 10 ÷ 20 W, where W is the width of the envelope of active channels. This minimum length corresponds to a number K of segments (K = Ls/d, with d length of a discretization segment) and their associated values of the considered variable;

◦

identify the discontinuity points, i.e., locations where the binary indicator (measuring the relevant attribute) changes from 0 to 1 or vice versa, and determine the distance D(-1) from the previous discontinuity (the notation D(-n) denotes the number of points separating previous discontinuity from the segment located n steps before current one);

◦

if the current point is not a discontinuity, then just keep the original value of the variable; if it is a discontinuity, instead, “look forward and backwards” to ascertain whether the discontinuity is due just to a local effect (“absence’, i.e., a “hole”, or “presence”) or it is the beginning of a new reach within the significant horizon, and assign the new value accordingly (see Figure 8).

→

All of this results in a sequence of value for the new binary variable (0,1) “presence” called “Holistic 1”.

Now, run the following procedure for the newly calculated binary variable Holistic 1:

◦

identify new discontinuities for this new variable;

◦

to each one, apply a forward moving averaging (on the K horizon) operator and, if the result is less than the threshold α, come back to a 0 value (absence), otherwise set it to 1;

◦

in this way, obtain a sequence of values for the new binary variable “presence” called “Holistic 2”;

◦

in this series, filter out all reaches characterized by absence, but shorter than the minimum significant length L_S;

◦

obtain thus a sequence of values for the new binary variable “presence” called “Holistic 3”;

◦

in this series, filter out all reaches characterized instead by presence, but shorter than the minimum significant length.

→

Result is the final sequence of values for the binary variable called “Holistic ok”.

This algorithm can be run for a “minimum significant length” L_S varying along the river, according to the fluctuating (generally increasing) width of the active channel envelope.

We are aware that this is just a first version that can certainly be improved (some steps are perhaps unnecessary); however, as shown in the next section, it performs quite nicely, at least in our case study.

3.2.3. Binary Case: Application to Magdalena River’s Islands

“Presence of islands” is an example of a binary attribute where a manipulation of the originally assessed values (1: presence, 0: absence of an island in each discretization segment) may be of interest if one wants to come up with an identification of reaches “with islands” as distinguished from those “without”.

From Figure 9, it can be observed that no reaches shorter than the established significant length L_S are left; but also that sometimes (particularly, in the first reach at the bottom of the first stretch at the left) isolated (tiny) islands appear in a stretch labelled as “no islands”; this is an expected result that is consistent with the way the algorithm is set up.

3.2.4. Reductionist→ Holistic Algorithm: The Categorical Case

A categorical variable assumes discrete values defined over a non-ordinal scale (i.e., it is not possible to rank them in any order). As an example, for the Cauca River (main tributary of the Magdalena) we created a variable denoted “current,” related to the type of surface flow (as recognizable from Google Earth imagery; this is a proxy for “character of the bed,” a variable that cannot be assessed directly for turbid rivers like these, where the bed is in general not visible). The variable “current” has seven possible categories (and corresponding discrete values on a non-ordinal scale):

0: rapids; 1: corrugated (long riffles that do not end in a pool); 2: pool and riffle; 3: streamed (stretch with visible lines of current and local turbulence, but softer than in previous cases); 4: rippled (softer than streamed); 5: plane surface (probably corresponding to planar bed: runs or glides); 10: not detectable (bad quality of images).

The philosophy for the holistic algorithm in this case is very simple: “do not leave windows shorter than L_S and assign to each stretch the prevailing value in it”. It can be noted that this algorithm is structurally “rougher” than in the binary case because there is not a clear criterion to choose which category to keep when a stretch is too short and hence we adopted a very simple and practical criterion that, although working acceptably well, could be well substituted by others.

Hence, the proposed algorithm works as follows:

Cycle I:

◦

Identify discontinuities (where switches from one category to any other occur).

◦

Where no discontinuity is present, the algorithm keeps the previously calculated value (as a result of this holistic algorithm; at the first step it assumes the original value).

◦

Where a discontinuity occurs:

-

“K steps forward prevailing”: the algorithm identifies that category which occurs with the highest frequency in the K steps downstream.

-

"Residual frontal window": analogous, but along just K–D(-1) steps forward (this is a forward moving horizon which progressively narrows down within a fixed window starting from the previous discontinuity; it is modified when a switch occurs to the next discontinuity). Its purpose is to consider what the prevailing value was previous to the current point.

-

When the prevailing value in the residual forward window coincides with that previous to the last discontinuity → current segment is just a “local hole” and as such the algorithm keeps the previous holistic value.

-

When the value is "not detectable" (one of the 7 possible categories), the algorithm just keeps it.

-

If not, the algorithm takes the same value of the prevailing category within the K steps forward.

Cycle II: removing residual segments shorter than K steps.

Here the following (arbitrary) criterion is adopted: “just make them uniform with the previous adjacent category”. Hence, the algorithm goes like this:

-

It identifies the discontinuities in the result from Cycle I.

Where no discontinuity is present, it keeps the previous value (the new one, produced by the algorithm; where a discontinuity occurs instead:

-

If the full distance from the previous discontinuity is smaller than K, the algorithm assigns the same value occurring before the discontinuity.

-

If not, it keeps the current value.

Again, the fitness criterion here is purely visual (Figure 10).

4. Conclusions

We can honestly say that the proposed systematic procedure and set of computer-aided tools (the ToolBox, some of whose tools we presented herein) do provide a significant benefit. An objective and systematic identification of reaches (addressed in Nardini et al., [18]), a reliable (although perhaps somehow cumbersome) confinement procedure, the reductionist-holistic tools, and the automatic, GIS based, grouping tool (not described in detail here: see Table A1 of Appendix A) provide all together a substantial support.

We are aware that the procedure reported in Table A1 and the specific algorithms and tools we developed (e.g., the confinement tool) still require a significant manual intervention and sometimes not straightforward exchanges between GIS and Excel; in other words, the proposed procedures are automatable, but not yet fully automated. All this, however, lies just at the operational level; what really matters is that the expert judgment is actually needed only in specific, well identifiable tasks: the definition of criteria like the Procedural Tree; the fine revision of the VB identification; the overall interpretation of the river behavior; and progressively less, the mapping of geomorphic units (but this is very likely to be soon performed by automatic procedures as discussed, for instance, in Piégay et al. [1]).

The conceptual clarifications and operational procedures, tools, and hints provided within the ToolBox and associated Manual can significantly support many practitioners—particularly in the Spanish-speaking world—and can help disseminate the adoption of the River Styles Framework, as well as the use of the Fluvial Corridor ToolBox.