An Integrated Approach to Riverbed Morphodynamic Modeling Using Remote Sensing Data

Abstract

River inlets, deltas, and estuaries represent delicate ecosystems highly susceptible to climate change impacts. While significant progress has been made in understanding the morphodynamics of these environments in recent decades, the development of models still requires thorough testing and data integration. In this context, remote sensing emerges as a potent tool, providing crucial data and the ability to monitor temporal changes. In this paper, an integrated approach combining remote sensing and morphodynamic modeling is proposed to assess river systems comprehensively. By utilizing multispectral or RGB optical imagery from satellites or UAVs for river classification and remotely derived bathymetry, echo sounder data for ground truth, and photogrammetric modeling of emerged areas, we outline a procedure to create an integrated and continuous digital terrain model (DTM) of a riverbed, paying particular attention to the wet–dry interface. This method enables us to identify the river centerline, its width, and its slope variations. Additionally, by applying a linear morphodynamic model that considers the spatial variability of river morphology commonly found in estuarine environments, it is possible to predict the wavelength and migration rate of sediment bars. This approach has been successfully applied to recreate the DTM and monitor the morphodynamics of the seaward reach of the Roya River (Italy).

1. Introduction

Small river basins frequently experience rapid flood events that can result in significant loss of life, destruction of resources, and damage to structures and infrastructure. These events have substantial economic, social, and psychological impacts. Over the past thirty years, the Liguria region in Italy has witnessed several intense flood events, including those in Sturla (1992), Varenna (1993), Verbone and Armea (2000), Chiaravagna (2010), Vara and Cinque Terre (2011), Bisagno-Fereggiano (2011 and 2014), Entella-Rupinaro (2014), Polcevera-Fegino and Stura-Orba (2019), and Roya (2020) [1].

Small and steep river catchments are highly vulnerable to extreme rainfall events, both due to their topography, which increases the probability of flash floods, and the lack of detailed hydrologic, hydraulic, and morphodynamic information needed to predict their stability under extreme floods and thus to manage these environments. Additionally, the rapid onset of floods in these areas complicates the design of protective and warning measures. A recent example in Liguria occurred in the Roya basin on the French–Italian border, which was hit by the Alex storm in 2020 [2,3], causing the Roya River to permanently alter its course and occupy an adjacent national road after an extraordinary flood event. This underscores the necessity of understanding hydraulic and morphodynamic processes in river basins to improve the protection and management of these fragile environments against flooding.

Over the past few decades, significant advancements have been made in understanding the morphodynamics of river courses and inlets [4,5]. However, developing and validating morphodynamic models requires data that accurately reflect the geometry of natural environments and the factors influencing their spatial and temporal evolution. While such data were once difficult to obtain, derived mainly from surveying cross-sections of the river plain by total station or GPS/GNSS satellite positioning, remote sensing techniques, such as photogrammetry or laser scanning from UAV or aircraft, now allow for the rapid and extensive acquisition of river morphology. In fact, these techniques enable the production of high-resolution, spatially continuous, and accurate data [6,7,8].

Even these remote sensing technologies face limitations. Photogrammetry cannot penetrate dense vegetation; it captures the height of the canopy rather than the ground beneath. On the other hand, Light Detection and Ranging (LiDAR) sensors, used by UAVs or aircraft, are effective in penetrating vegetation canopies [9]. However, given the high cost of LiDAR sensors, photogrammetry remains a valuable tool for reconstructing the emerged areas of a river DTM.

Submerged areas are traditionally surveyed using Single-Beam or Multi-Beam Echo Sounders (SBES or MBES) mounted on boats or uncrewed watercraft. They produce highly accurate bathymetric data [10,11] but require a significant investment in equipment, which must be carefully calibrated, and considerable time for data collection and processing [12,13]. In the case of shallow waters, bathymetry can also be estimated by photogrammetry or bathymetric LiDAR after appropriate correction for water refraction [14]. In addition, multispectral optical images from satellites, aircraft, and UAVs have also recently been used to obtain remotely derived bathymetry (RDB) products [15,16,17,18]. However, in the case of unclear waters, estimating bathymetry from remote sensing can prove technically challenging [19,20].

Several methods have been proposed to obtain digital terrain models (DTMs) of emerged and submerged areas in coastal or riverine environments by integrating remote sensing data from various sources. In fact, a complete description of such environments is essential for many applications [21,22,23,24]. Integration is already an issue within the individual data acquisition domains, without necessarily referring to the further integration between bathymetry and topography, as addressed by Lewicka et al. [25]. A first method aiming to create a continuous surface for bathymetry and topography in coastal areas was proposed by Gesch and Wilson [26]. While the premises of their study share some similarities with ours, their focus on a shallow coastal environment distinguishes it from our riverine context. They acquired high-quality data and successfully interpolated along the coast. Integrating data in fluvial environments may require additional attentions, as fluvial topographic and bathymetric data are often collected using very different instruments. A comprehensive DTM creation process was proposed by Schäppi et al. [27] relying on riverbed cross-sections, which were not available in our dataset. Note that the accuracy of cross-sections is depth-dependent, since they are usually performed by walking in the riverbed. While the method is suitable for fluvial environments, it requires specific adaptations to handle complex morphologies and significant topographic changes. Additionally, the quality of the model depends on the density of cross-sections and the accurate definition of breaklines, especially in dynamic areas such as river bends or reaches subject to frequent morphological changes. Karaki et al. [28] explored more experimental alternatives, integrating a UAV carrying an aerial optical sensor for photogrammetry, with an ASV (Autonomous Surface Vehicle) equipped with underwater acoustic sensors. That solution required specialized instrumentation, which usually raises costs. In this regard, our goal was to develop a highly adaptable method for fluvial environments that could effectively integrate heterogeneous datasets, even those not initially intended for merging, while maximizing automation and replicability.

This paper shows the research performed so far within the European-funded RAISE1. It aims at integrating morphodynamic modeling with digital modeling derived from image analysis, for enhanced river management. Images enable continuous monitoring of the river, while the morphodynamic model facilitates management by offering quick and user-friendly simulations to address critical questions, such as the following ones: Will dredging in one area lead to deposition elsewhere? Will bedform migration cause shifts in the excavation site, and if so, how quickly? What long-term changes can be anticipated if river flow increases?

The contributions of this study are as follows: First, we propose a method for generating a comprehensive DTM of the entire riverbed by integrating image analysis, taking care of its automation and replicability for future application to different case studies. This is crucial for monitoring riverbed changes over time and providing valuable information that can be used to test morphodynamic models. Second, we use the morphometric data derived from the DTM, such as channel width, bar length, localized scours, and deposition areas, to calibrate a simple morphodynamic model and assess the accuracy of its predictions.

Concerning the DTM generation, the submerged morphology is computed using a Remote-Derived Bathymetry (RDB) algorithm [17,29] applied to multispectral optical imagery, while the emerged sections are measured using UAV photogrammetry. Concerning morphodynamic modeling, we refer to the model proposed by Ragno et al. [30], based on a linear stability analysis. It offers the advantage of predicting, at a low computational cost, whether bars may form and if so, the type of bars (alternate or central) that typically emerge in the early stages of development, along with their wavelength and migration rate. Moreover, the model represents the forefront of research in the literature on linear stability analyses of river bedforms, as it accounts for the effects of seaward boundaries and width variation.

This integrated approach is applied to create a DTM and monitor morphodynamics in the seaward stretch of the Roya River (Italy). The present findings indicate a qualitative alignment between the theoretical model and remote sensing data, encouraging the further development of the present integrated approach to improve understanding of river morphology, hydrodynamics, and sediment dynamics, to support informed decision-making for sustainable river management. Additionally, the proposed method emphasizes replicability and cost-effectiveness, contributing to ongoing river monitoring and enhancing datasets for model refinement.

The structure of the article is thus organized as follows. Section 2 offers an outline of the state of play of the principles and methodologies behind morphodynamic modeling, along with insights into the creation of the DTM. Section 3 outlines the materials, including the study area and data collection information. Section 4 explains the methodology used for DTM generation and morphodynamic modeling. Section 5 presents the results, followed by a discussion of the study’s advantages and limitations in Section 6. Finally, Section 7 provides the conclusion and future work.

2. Background

2.1. Morphodynamic Modeling

From a morphodynamic standpoint, rivers may conceptually be modeled as systems in a state of quasi-equilibrium, i.e., such that flow perturbations relative to some equilibrium state let the system evolve on spatial-temporal scales much larger than the hydrodynamic scale. This quasi-equilibrium pertains to the cross-sectionally averaged bed elevation rather than its local value. Indeed, in the real world, rivers display spatial variations and temporal fluctuations in the flow properties acting on various scales, like free and forced bedforms, which arise from bed interface instabilities and channel geometry deviations, respectively. Moreover, larger scales involve variations in water discharge and sediment flux due to events like floods and seasonal changes. Despite these perturbations, the concept of morphodynamic equilibrium, viewed in an averaged sense, remains insightful [5].

Traditionally, equilibrium in geomorphology is linked to the notion of formative discharge [31], which essentially assumes that the unsteady forcing on the river is morphologically equivalent to some effective steady forcing. The initial goal of morphodynamics was thus to define this equilibrium state in a simple configuration: a cohesionless straight channel with a constant slope under steady flow and sediment supply. Despite its existence, mechanistic approaches have revealed that the equilibrium configuration is not stable. Notably, bed perturbations lead to the formation of large-scale fluvial sedimentary patterns known as bars.

The so-called “free bars”, resulting from spontaneous bed state instability, exhibit large-scale deposition bumps and scour holes on alternate channel sides and are typically migrating features. These bars impact navigability, enhance bank erosion, and interact with engineering structures. For these reasons, over the past 50 years, extensive literature has emerged on river bars, including theoretical studies [32], experimental investigations [33], field observations [34], and numerical modeling [35]. These studies show that bars form spontaneously in both gravel and sand bed rivers when the channel’s width-to-depth ratio exceeds a critical value. Bars are downstream-migrating features, often appearing as single rows (alternate bars) or multiple rows in wider, shallower streams. A linear analysis predicts bar wavelength and migration speed, while a nonlinear analysis is needed to determine their amplitudes [36].

Compared to the pure riverine case, the morphodynamics of estuaries, in which both the fluvial and tidal forcing concur to determine the flow field, has received comparatively less attention. Most research has concentrated on understanding the long-term equilibrium profile [37,38], with some recent numerical simulations exploring the morphodynamic evolution of large-scale estuarine systems [39]. However, a complete study dealing with bars’ formation, evolution, and equilibrium configuration in these peculiar transitional systems is still absent. For this reason, here, we refer to the model of Ragno et al. [30], who faced the problem of the embryonic formation of free bars in weakly converging tidal channels and estuaries mathematically within a unified two-dimensional shallow water and sediment mass balance model. The main physical parameters employed in the model are recalled in Section 4.2 for its applications.

2.2. Digital Terrain Model of a River

Morphodynamic modeling frequently necessitates data describing the river environment from both morphometric and hydrodynamic perspectives. An accurate DTM representation of the riverbed, with continuity between the emerged and submerged parts, is necessary both as input data and for model validation [40]. It allows one to extract morphometric information, such as channel width, bar length, localized scours, and depositions. Moreover, information on sediments in the riverbed and the possible presence of vegetation are fundamental for the definition of the roughness coefficient, which influences water velocity and discharge. Obtaining such data is usually challenging and requires costly measurement campaigns. Usually, descriptive elements of the riverbed are sourced from ground-based techniques, such as total stations or GPS/GNSS positioning, involving traversing the river plain and requiring each point to be approached individually [41]. This time-consuming methodology requires the watercourses to be accessible and the water to be shallow. Furthermore, the acquired data are discrete. Recent advancements in technology have significantly streamlined the processes of data acquisition. Nowadays, various remote sensing techniques enable the production of high-resolution, spatially continuous and accurate data, following meticulous calibration, processing, and analysis [6,7,8].

Tools valuable for remote sensing studies of river morphology include aerial and satellite imagery [42]. Satellite images offer a synoptic and recurrent viewpoint on the area of interest, providing automatic acquisition [43]. On the other hand, aerial imagery can be produced by UAVs, yielding high spatial resolution, flexibility, capacity to gather data at different times of the day, and relatively lower cost compared to traditional platforms [44]. UAV photogrammetry enables extensive and rapid surveys, offering valuable information for modeling the non-vegetated emerged sections of the riverbed. However, in areas with dense vegetation, it captures the height of the canopy rather than the underlying terrain.

Light Detection and Ranging (LiDAR) is a remote sensing technology used from UAVs or aircraft that measures distances by emitting laser light and analyzing the reflected signal. It is effective for surveying both dry and wet areas. One of LiDAR’s key advantages is its ability to penetrate vegetation canopies [9], providing detailed information about ground-level data even in densely vegetated areas. When LiDAR is used for topographic surveys, it typically operates at near-infrared (NIR) wavelengths (1064 nm) [45,46] and the signal cannot penetrate water, thus providing a direct reference for water surface elevation [47]. In contrast, bathymetric LiDAR uses green light that travels through the water column up until reaching the submerged riverbed. However, this technology can face challenges related to transient water turbidity, which affects data quality in shallow water environments. Separating surface reflections, water column effect, and bottom reflections in LiDAR-based surveys can also prove technically complex [20]. Given the high cost of LiDAR sensors, photogrammetry is a technique most widely used for reconstructing the emerged areas of a river DTM.

Submerged areas are traditionally surveyed using acoustic sonars mounted on boats or uncrewed watercraft. Single-Beam Echo Sounders (SBES) and Multi-Beam Echo Sounders (MBES) are well known for their ability to produce highly accurate single-point or continuous bathymetric data [10,11]. However, achieving this level of precision requires significant equipment investment, careful installed device calibration, and considerable time for surveying and data processing [12,13]. Additionally, using these technologies in shallow waters can be problematic due to breaking waves and submerged obstacles, limiting their overall applicability [48].

In the case of very shallow and clear water, bathymetry can be estimated from photogrammetry after appropriate correction for water refraction [14]. In the case of unclear water, an appropriate relationship between turbidity and brightness factors can be found [19].

In recent decades, optical multispectral imagery from remote sensing platforms, such as satellites, aircraft, and UAVs, has also been used to obtain Remote-Derived Bathymetry (RDB) products [15,16,17,18]. These techniques rely on analyzing the interactions between light at various wavelengths and the water column, which are influenced by water constituents and bottom materials [49]. There are two primary approaches for estimating water depth based on optical band attenuation in RDB [15]. Statistical approaches entail correlating the reflectance of one or more suitable spectral bands with bathymetry, without considering the physics of light propagation in water or environmental factors [49]. Their main advantage is their ease of application. Conversely, physics-based approaches incorporate the dynamics of light propagation in water. These models can be solved using semi-empirical and analytical methods. Semi-empirical methods are founded on specific assumptions about radiative transfer and light attenuation in water and are calibrated with in situ data [15]. Notable semi-empirical methods include those developed by Lyzenga [50] and Stumpf [29], which have been applied directly or adapted for deriving bathymetry in various studies focusing on fluvial environments [14,51,52,53]. Analytical methods depend on understanding the biophysical and optical properties of water constituents. While analytical methods can achieve high accuracy, they are complex to implement and require numerous in situ parameters and precise atmospheric correction inputs [49].

3. Materials

3.1. Study Area

Roya is a river in the French–Italian Maritime Alps (43°48′21″ N, 7°35′33″ E), flowing into the Ligurian Sea in Ventimiglia (Figure 1).

It is 59 km long and has an area of drainage basin of about 660 km². The source of the Roya River is located in Colle di Tenda, at 1908 m AMSL, within the French municipality of Tenda. It exhibits typical characteristics of an alpine river in its upper course, with a steep and rugged bed traversing narrow gorges. As it descends towards the valley, it flows through several localities, including Tenda, San Dalmazzo di Tenda, and Fontan, and is joined by tributaries such as the Bevera stream. The hydrology of the Roya River is characterized by a torrential regime, with substantial flow rates throughout the summer due to its alpine watershed. Nevertheless, the river is prone to sudden and intense flooding, particularly during autumn and spring. A notable event occurred in October 2020, when heavy rainfall caused by the Alex storm [2,3], together with the release of water from the Casterino dam, resulted in severe flooding and widespread damage along the river course, significantly affecting Ventimiglia, where the river mouth is located. It can therefore be inferred that the river was chosen for this study due to its societal relevance. The reach under investigation begins at the river mouth and extends up to the confluence with the Bevera Stream, spanning a linear distance of 4.67 km in a north–west direction, as illustrated in the red rectangle of Figure 1. This reach was selected as it offered the best balance between the spatial coverage of the available data, the gradual variation in channel dimensions moving upstream from the mouth, and the requirement for a sufficiently large area to effectively evaluate the resulting digital model.

3.2. Data Collection

Datasets employed in this study are described in

Table 1. Datasets.

Data Type	Source	Acquisition Date	Features
Aerial orthophoto	AGEA	13 August 2019	0.2 m spatial resolution
Photogrammetric digital model	Gter	Summer 2019	0.4 m spatial resolution
SBES points	CNR-INM/Gter	5 November 2019	42,585 sampled points covering an area of approximately 0.04 km2

.

• The aerial orthophoto (Figure 2) available on the Liguria Region Geoportal (Italy)2 was acquired in 2019 by RTI CGR Spa/e-Geos S.p.A under a commission received from AGEA. The image comprised 4 spectral bands (red, green, blue, NIR) with a spatial resolution of 0.2 m. The image was used to classify the riverbed and calculate the bathymetry.
• The photogrammetric survey (Figure 3), performed by Gter Ltd. during the summer of 2019, produced a digital model with spatial resolution of 0.4 m. In this work, the photogrammetric digital model was used to extract bare-soil areas’ elevation and to create a surface that allowed bathymetry to be converted from depth to elevation.
• The SBES survey, performed within the seaward reach of the Roya River (Figure 2), relied on an Autonomous Surface Vehicle (ASV), referred to as SWAMP. It was equipped with a Microstrain 3DM-GX3-35 GPS + AHRS (attitude heading reference system), and an Echologger ECS400 single-beam Sonar (with a working frequency of 200 kHz) for river-bottom range measurement [54]. The survey, carried out on 5 November 2019, sampled more than 42,000 points over an area of approximately 0.04 km². In this study, SBES points were used to calibrate and validate the RDB model.

4. Methods

4.1. DTM Computing Method

The UML (Unified Modeling Language) activity diagram of Figure 4 illustrates the workflow to reconstruct the DTM of a riverbed from remote sensing data, described in detail in the following subsections. The initial step involves identifying whether sections of the riverbed are above or below the water level (Section 4.1.1). Emerged areas are then categorized by riverbed type (i.e., bare soil or vegetation), while an algorithm is applied to submerged ones to retrieve bathymetry (Section 4.1.2). Subsequently, Elevation Point Clouds (EPCs) are generated for each section, based on the specific data sources available. Lastly, EPCs are interpolated to reconstruct a continuous DTM of the riverbed (Section 4.1.3).

4.1.1. Image Classification

Classifying the elements within the riverbed involved processing the multispectral image through a series of sequential steps, providing the basis for all subsequent workflow activities.

First, the internal portion of the riverbed bounded by the lateral banks had to be identified. To this end, we referred to a land use map from the Ligurian Geoportal, which was manually refined to single out the area of interest and mask certain anthropogenic features, such as overpasses, bridges, and side walls, as shown in Figure 2.

Subsequently, distinguishing emerged and submerged sections of the riverbed was accomplished by analyzing the near-infrared band. Since water has a distinct spectral signature in this wavelength range compared to bare soil and vegetation [55], by applying an appropriate threshold value to this band, it was possible to differentiate wet and dry surfaces.

Through a supervised process implemented in a GIS environment, the emerged sections were further classified between vegetated areas and those comprising bare soil, which in the seaward reach of the Roya River mainly included small pebbles. Training maps representative of the two different homogeneous classes within the image were manually defined, and the corresponding spectral signatures were computed. These served as input to a maximum likelihood classifier algorithm3, which assigned each pixel previously marked as dry land to the class to which it had the highest probability to belong.

Similarly, a supervised classification was also employed to identify shadows cast on the water surface by vegetation and crossing structures. This was necessary because shadows have a negative impact on the calibration of the model used to calculate bathymetry, as illustrated below.

4.1.2. Remote-Derived Bathymetry

Bathymetric data were derived from optical multispectral imagery extending the semi-empirical method proposed by Stumpf et al. [29]. That method is based on the principle that light intensity decreases exponentially as it passes through the water column due to absorption, with the rate of attenuation depending on the wavelength. It relies on a set of in situ bathymetric observations to correlate water depth with the ratio of log-transformed reflectance values from appropriate spectral bands. In this study, drawing on suggestions from other authors who analyzed variations to Stumpf’s linear solution, offering better adaptability and accuracy in capturing the observed data’s behavior [56,57], we extended Stumpf’s approach by adopting a three-parameter exponential model:

$$Z=\mathit{exp}\left\{\mathit{a}{\displaystyle \frac{\mathit{ln}\left[\mathit{n}{\mathit{R}}_{\mathit{w}}\left({\lambda }_i\right)\right]}{\mathit{ln}\left[\mathit{n}{\mathit{R}}_{\mathit{w}}\left({\lambda }_j\right)\right]}}+\mathit{b}\right\}+\mathit{c}$$

(1)

where Z is the water depth, $R_w$(${\lambda }_i$) and $R_w$(${\lambda }_j$) are the reflectance values measured across the spectral bands having ${\lambda }_i$ and ${\lambda }_j$ as their central wavelengths, a, b, and c are calibration parameters, and n is a fixed constant ensuring the ratio always remains positive.

The high-density SBES data were resampled to align with the spatial resolution of the orthophoto, preventing multiple depth values from being assigned to the same pixel. This operation reduced our initial 42,585 ground-truth points to 9752. Finally, resampled points were randomly divided into two separate datasets: one for calibration and one for validation of the RDB model, each comprising half of the available points.

Model calibration was performed following the nonlinear least squares method. An initial set of parameters that provided the best-fit regression between the calibration dataset and the image band ratios was computed. Both the blue/red and blue/green band ratios were examined, with the spectral bands masked to show only the submerged sections of the riverbed. The blue/red ratio was ultimately selected due to its superior performance in shallow waters. The coefficient of determination (R²) was computed as a proxy for the quality of the model calibration. This statistical measure indicates the percentage of variance in the dependent variable that can be predicted by the independent variable. A higher R² value indicates a better fit, meaning the model more accurately captures data trends. Outliers were removed according to Chebyshev’s theorem [58], also known as the three-sigma criterion, which excludes values more than three standard deviations from the sample mean. The model was then recalibrated, and a validation procedure was carried out to evaluate its robustness. This included comparing the estimated depths with the actual validation dataset and calculating a series of statistical parameters: Root-Mean-Square Error (RMSE), Mean Absolute Error (MAE), Bias Average (BIAS AV), and Bias Standard Deviation (BIAS STD). The RMSE quantifies the average squared deviation between predicted and observed values, with lower RMSE values signifying higher accuracy of the model. The MAE captures the average magnitude of errors in predictions ignoring their direction. The BIAS AV refers to the average residuals, computed as the difference between control depths and the estimated RDB, showing whether the model overestimates or underestimates the target variable. The BIAS STD indicates the variability of these residuals around the BIAS AV, with lower values indicating more consistent predictions.

4.1.3. DTM Reconstruction

EPCs were generated for each riverbed type (i.e., bare soil, water, and vegetation) based on the specific data sources available. Then, EPCs were interpolated and eventually smoothed to obtain a discontinuity-free DTM that encompassed both emerged and submerged areas.

Identification of EPCs for submerged section. While geomatics instruments provide direct elevation readings for emerged areas, submerged areas, which are surveyed using echo sounders or RDB, are commonly described by depth values rather than elevation measurements [59]. Depth is measured as the vertical distance between the free-water surface and the river bottom. Thus, by knowing both the depth values and the elevation of the water surface, which follows the gradient of the river from the mouth to the source, the river bottom elevation can be properly derived.

The proposed method calculated water depth from multispectral images and estimated the water surface elevation using bare soil points located just next to the water as indirect measurements. This process began by identifying the river’s centerline, which could be carried out either manually through digitization or by using Voronoi polygons [60] if a set of points along the river-banks was available. Once the centerline was established, a series of points was spaced along it at 1-m intervals. Next, a buffer zone 0.5 m wide was created around the water areas. Within that buffer, points classified as bare soil were identified, as these points were considered representative of the water surface height. The elevation of each point along the river’s centerline was then assigned based on the closest bare soil point within the buffer zone. Finally, using Voronoi polygons, the elevation values assigned to the centerline were extended across the entire riverbed, allowing for a continuous representation of the water surface elevation based on the surrounding terrain. This approach leveraged the spatial relationship between bare soil and water to generate a detailed model of the river’s surface.

This method created a stepped surface, which was then filtered using a SAGA Simple Filter tool [61] to produce a smoother result mirroring the height increase from the river mouth to its source. The riverbed elevation was computed by subtracting depth values from the water surface height.

Identification of EPCs for the emerged section. Elevation values for bare soil were directly extracted from the photogrammetric point cloud. However, since photogrammetry does not provide detailed elevation information for the ground beneath the canopy, the proposed approach used interpolation techniques to estimate the EPC beneath the canopy. The process began by creating a buffer zone (0.5 m) around the vegetated areas. By subtracting the original vegetation area from that buffer, a thin strip outside the vegetation was obtained. This strip included water points, with elevations calculated as the difference between the water surface and the RDB, and bare soil points derived from photogrammetry. The elevations of water and bare soil points were then interpolated using a Triangulated Irregular Network (TIN) to produce a surface that approximated the ground elevation beneath the vegetation. Finally, that interpolated surface was converted from raster to vector format to generate an EPC in the vegetated areas.

Generation of DTM from all the EPCs. EPCs for each class (i.e., bare soil, water, and vegetation) were consolidated into a single dataset. TIN interpolation was then applied to generate a continuous surface. That method is particularly effective for filling in large no-data gaps, such as those beneath bridges and walkways. To further smooth the resulting surface, it was converted from a raster grid to a vector point cloud, and cubic spline interpolation was applied. If necessary, an additional smoothing step may be performed using a Gaussian filter [61]. The final result of this process was a continuous, smooth surface that provided elevation values for both the emerged and submerged sections of the river.

4.2. Morphodynamic Method

In this application, the morphodynamic model for bar formation in micro-tidal estuaries and tidal channels developed by Ragno et al. [30] was utilized. The main features of the model are hereafter summarized, while interested readers may refer to [30] for a detailed formulation.

The model operates within a two-dimensional framework, using a depth-averaged shallow water formulation commonly employed in many bar models. The flow is assumed to occur in a well-mixed tidal stream, over a cohesionless bed composed of homogeneous sediment characterized by a median grain size $d_s$, and a tidal sea as the seaward boundary. The effect of wind waves is neglected, and the banks are assumed to be non-erodible. A Cartesian reference frame is used, with x and y indicating the longitudinal and transverse coordinates, respectively, with $x=0$ at the tidal sea and extending landward. Spatial variations in channel geometry are considered, as depicted in Figure 5.

In natural tidal environments, channel width typically increases gradually toward the mouth, which can be approximated by an exponential law:

$$B=B_u+(B_m-B_u)exp\left(-{\displaystyle \frac{x}{L_c}}\right)$$

(2)

where $L_c$ is the convergence length, $B_m$ is the half-width at the mouth, and $B_u$ is the upstream half-width not affected by tidal processes.

The model relies on several spatial scales: the planimetric dimension of bars related to channel width, the convergence length governing the rate of smooth widening toward the mouth, and the backwater length $L_b$, defined as the ratio between uniform flow depth and bed slope, crucial for studying coastal processes like stratigraphic patterns and avulsion phenomena. Analyses have shown that B is typically much smaller than both $L_c$ and $L_b$, which are of similar magnitude. This identifies two spatial scales for studying bedform formation: the short bar scale and the long channel scale. The governing equations are the St. Venant equations for shallow water flow, coupled with the Exner equation for sediment mass conservation. Closure relations are needed to express shear stress and sediment transport rate in terms of the governing variables. This is achieved through a “local equilibrium” approximation, justified by the slowly varying flow field. Shear stress and sediment transport rate are computed based on local flow variables, with the friction coefficient estimated using the Manning–Strickler formula and sediment transport rate determined using the Engelund et al. [62] formula.

To solve the system, proper boundary conditions are imposed: the channel walls are impermeable to both flow and sediment fluxes, the tidal sea induces a small oscillation of the free surface at the mouth, and normal flow conditions with steady freshwater discharge per unit width and bed slope are assumed upstream.

According to Ragno et al. [30], bar formation in tidal streams is investigated using a linear stability analysis, a mathematical tool to examine conditions where the system’s reference state loses stability to small perturbations. For an infinitely long straight channel with constant upstream water discharge, the equilibrium state is a uniform turbulent flow over a sloping bed. In tidal streams, channel widening and tidal motion complicate this basic flow, making it a function of time and space [38]. The linear analysis examines the stability of the basic state to small perturbations in flow and sediment transport. A linear “normal mode analysis” is performed, expanding perturbations in the Fourier series and investigating each mode independently. The linearized differential system is solved numerically using a Runge–Kutta scheme, starting with a small random perturbation of the bed topography.

The morphodynamic model requested input morphometric information, such as river width and slope, which were extracted from the riverbed DTM following the workflow outlined in Figure 6. In particular, the river centerline was first derived from the DTM. This allowed us to identify the curvilinear longitudinal axis of the river and each transversal cross-section of the river.

Each cross-section provided a local value of the river width (as the distance between the two cross-section banks) and of the bottom elevation (by averaging transversally along the cross-section the bottom elevations).

Once the distributions of the river width and cross-sectionally averaged bottom elevation along the curvilinear longitudinal axis were known, these data could be analytically fitted to extrapolate key parameters for the morphodynamic model, such as the mean river slope, inlet width, uniform flow width, and convergent length.

Additional field data directly extrapolated from the DTM included bar patterns, which were identified by rectifying the bottom profile and cutting back the cross-sectionally averaged bed profile. These data were used to validate the morphodynamic model by comparing them with the corresponding information generated by the model.

5. Results

The methods described above were applied to create the DTM and study the morphodynamics of the seaward reach of the Roya River, from its mouth to the confluence with the Bevera Stream. To present clear and easily interpretable figures, the DTM process results are only shown for the estuarine terminal portion of the river. Conversely, the results of the morphodynamic model cover the entire analyzed reach of the river.

5.1. DTM Computing

The results of the Roya riverbed classification, obtained by processing the multispectral orthophoto in GRASS GIS, are displayed in Figure 7a. Submerged sections of the riverbed are shown in blue, vegetated areas in green, bare soil composed of pebble-sized sediments in yellow, and shadows cast on the water surface in red.

Figure 7b presents the bathymetric map generated using the band ratio method. The robustness of the RDB model was confirmed by the calibration and validation statistics in Figure 8. The high R² value of 0.868 suggested a strong correlation between depth values and the band ratio. In addition, the low RMSE value of 21 cm denoted that the model achieved a good level of accuracy in depth estimation, effectively capturing the variations in the riverbed topography.

Figure 9 presents the results of the multi-step process used to create a stepped surface, which enabled the transition from a depth-based reference to an elevation-based one for the submerged sections of the riverbed. Figure 9a shows the river centerline, identified by constructing Voronoi polygons from the points along the river-banks. The pink dots represent bare-soil photogrammetry points located within the buffer zone surrounding water areas, whose elevation was assigned to the nearest points along the river centerline (yellow dots). Voronoi polygons were constructed using both the centerline and banks points. Figure 9b shows the final water surface outcome. It is important to note that the elevation does not increase linearly from downstream to upstream, particularly near the river mouth, where significant deposition occurs. In that region, the high variability in the elevation of bare soil points close to the water led to uncertainty in determining the water surface elevation. Moreover, since the riverbed elevation was calculated by subtracting RDB from the water surface elevation; where these two quantities are comparable in magnitude and uncertainty, as in the mouth area, this may compromise the ability of the proposed method to correctly reconstruct the elevation of the submerged area. Differently, moving upstream, the reconstruction of the submerged areas of the riverbed was more accurate because the elevation of the water surface increased and, at the same time, the uncertainty in its definition decreased due to the greater homogeneity in the levels of the bare soil points close to the water.

In Figure 10, the results for bare soil areas and for interpolation in the vegetated areas are shown. Figure 10a illustrates the elevation data for bare soil areas. This result was strictly dependent on the photogrammetric survey’s characteristics since the values were imported from the digital surface model without any processing. In Figure 10b, points within the buffer created around the vegetation are visible. These points, categorized as either water (in light blue) or bare soil (in orange), had their elevation values determined by subtracting the RDB from the water surface or by using photogrammetry, respectively. They were interpolated to reconstruct the terrain elevation beneath the canopies, as shown in Figure 10c. Here, some peaks in the digital model are present, mainly due to the high elevation in bare soil points at the boundary of the river plane. This outcome was influenced by challenges in the classification process: isolated points may have erroneously been identified as bare soil, resulting in their elevation being used in the DTM construction, when in fact they represented the top of the vegetation canopy. Consequently, the final result hinged on the photogrammetric survey’s ability to accurately capture the elevation of bare soil near steep banks or dense vegetation, as well as the precision of the initial classification of different riverbed types.

Figure 11 illustrates the final DTM, in both 2D (Figure 11a) and 3D visualization (Figure 11b). Since it was created through the interpolation of EPCs, the resulting surface was continuous, even in areas beneath bridges and overpasses, indicated by black-striped zones in Figure 11a, where elevation data were generally missing. This continuity made the DTM well suited for extracting the necessary information for morphodynamic modeling.

5.2. Morphodynamic Results

The 1D morphodynamic equilibrium of the terminal reach of the Roya River was analyzed through the DTM data obtained with the procedure described above. Moreover, other input data of the morphodynamic model were sediment characteristics and typical values of the flow rate, detailed by Colombini and Bolla Pittaluga [63].

The basic state (i.e., the cross-sectionally averaged bottom elevation) and the width distribution of the seaward reach of the river were derived by identifying the longitudinal axis of the river from the continuous DTM.

Figure 12 presents both the elevation of the mean bottom profile (represented by the green dots) and a constant slope profile (dashed line) that best fit the observed data. The real data points suggested a relatively uniform trend in the channel’s elevation, allowing for a simplified interpretation of the riverbed’s morphology. As illustrated, the channel’s slope could be approximated as constant, with an estimated value of 0.0052. This approximation of a constant slope allows for the application of theoretical models which often rely on idealized geometries. Nevertheless, as Figure 12 shows, it is important to acknowledge that real river systems often exhibit local fluctuations with respect to the mean slope due to sediment deposition, erosion, and natural or human-induced modifications to the channel. While these small-scale deviations are not captured by the constant-slope assumption, they may be incorporated into more detailed models or addressed in future studies where local morphodynamics are of particular interest. For the purposes of this analysis, however, the constant slope provided a reasonable first approximation that allowed for the study of the main river morphodynamic processes.

Similarly to Figure 12, Figure 13 shows the distribution of the river width along its longitudinal axis. As commonly observed in natural environments, the river width was characterized by a slight increase toward the mouth in its seaward reach (within the first km from the river mouth). Nevertheless, Figure 13 shows that the river width undergoes several fluctuations proceeding landward, on average slightly increasing its width. This peculiar behaviour may be associated with the fact that most river-banks are not natural but artificially fixed. For this reason, as a first approximation in this simulation, we considered the river width B constant and equal to an average value of 144 m (infinite convergence length $L_c=\infty$), leaving possible investigations on the effect of width variations to the future.

Additional input parameters for the model were the formative discharge, Q = 475 ${\mathrm{m}}^3$/s, the mean sediment diameter, $d_s$ = 29 mm, and the mean Strickler parameter of 30 ${\mathrm{m}}^{1/3}$/s.

Figure 14 illustrates the temporal evolution of the maximum bed perturbation amplitude for the first two lateral modes within the terminal widened reach, which spanned approximately the last 5000 m. The results indicate that linear theory identified the first mode (i.e., alternate bars) as the most rapidly growing mode. This finding aligns with the bedform pattern observed in Figure 15.

Figure 15 is also evidence that the considered river reach was not entirely straight but rather included bends along its course. However, the seaward half of the river, spanning from meter 0 to around meter 3600, was predominantly straight. This section of the river displayed minimal curvature, allowing it to be treated almost as a straight channel in preliminary analyses. The bends were confined to the landward half of the river reach, specifically between meter 3600 and meter 4700, covering an entire meander wavelength. However, recognizing that bends in the landward portion affected sediment transport, flow velocity distributions, and the formation of riverbed features like forced bars, the relatively straight seaward reach supported, as a first approximation, the employment of theoretical models referring to straight channels.

To better compare the real and the predicted bar wavelength, Figure 16 shows the longitudinal profiles of the bottom along the banks. The resulting bedform pattern obtained by subtracting the mean bottom elevation from the real one is reported in Figure 16a, suggesting a clear alternation of pools and riffles along the rivers, as simulated by the model (Figure 16b).

Field data (Figure 16a) indicated that bar wavelengths ranged between 700 and 900 m within the two kilometers from the mouth, decreasing to nearly half that distance in the third kilometer. Beyond 3 km from the inlet, no clear bar wavelength was detected.

The mathematical model (Figure 16b) provided similar results, with bar wavelength ranging about 800–900 m in the first kilometer from the mouth, and ranging about 600 m in the following two kilometers.

Some lag between the theoretical and actual bar profiles was expected, as bars are migrating features, and the model may not precisely capture their exact positions. Moreover, the theoretical model assumed a constant width, so some variation in the bar wavelengths was expected from that assumption. However, there was a general agreement between the mathematical model and the field data.

The two-dimensional plots of the corresponding bed topography are reported in Figure 17. In particular, Figure 17a displays field data directly extrapolated from the DTM, after rectifying the bottom profile and cutting back the cross-sectionally averaged bed profile according to the procedure described in Section 4.2. The typical rhythmic sequences of regions of scour (pools) and deposit (riffles) of alternate bars can be easily detected. It is important to note that the bed interfacial waves had wavelengths scaling with channel width and amplitudes scaling with flow depth. A similar alternate pattern of pools (red regions) and scours (blue regions) was predicted by the analytical model (Figure 17b). Although the linear model could not predict bar amplitudes, the comparison was still satisfactory, as it accurately forecast a bottom pattern with wavelengths closely matching the real ones.

6. Discussion

6.1. RDB Exponential Model

The RDB model used in this study to reconstruct the bathymetry of submerged river sections differed from the one proposed by Stumpf in the regression law applied. In the context of our study, the exponential solution showed greater flexibility in capturing data trends compared to Stumpf’s linear model, as indicated by the R² values in

Table 2. Statistical values resulting from the calibration and validation of the linear and exponential RDB band ratio models.

Model	R2	RMSE (m)	MAE (m)
Linear	0.840	0.233	0.172
Exponential	0.868	0.211	0.157

, which showed an improvement of about 3.3%. This choice positively impacted model validation as well. The RMSE and MAE values decreased by approximately 10%, demonstrating that the exponential model offered greater accuracy.

6.2. DTM Transect Analysis

Two transects of the riverbed were analyzed to evaluate the continuity and reasonableness of the final surface by comparing it with the geometries of the available reference data, as shown in Figure 18. The top section of the graphs indicates the riverbed classes (water, vegetation, and bare soil) along the transects. Corresponding reference lines are displayed for each segment: SBES data (light blue) serve as a reference for water sections, while the photogrammetric digital model (purple) provides a reference for vegetated and bare soil areas. The red line represents the final DTM elevation along the transect. In both transects, the DTM profile closely matched the reference data in the water and bare soil regions. In vegetated areas, discrepancies arose due to photogrammetry capturing the top of the canopy, whereas the DTM reconstructed the ground level through interpolation. Nevertheless, the interpolated surface still provided a realistic and accurate representation of the riverbed.

6.3. DTM Reliability Map

The DTM can be accompanied by a reliability map, which visually communicates the estimated confidence in the final data by highlighting the characteristics of the different data acquisition methods used for each riverbed section. Photogrammetric surveying applied to bare soil areas can be expected to yield centimeter-level precision [64], while RDB for water regions typically achieves decimeter-level accuracy [65]. In our study, RDB yielded 21 cm accuracy (RMSE value). However, the final elevation value for the water class was given by subtracting the RDB from the water surface height, whose accuracy could be lower than the decimetric level. In interpolated areas, such as those affected by vegetation, bridges, or shadows, the reliability of the DTM varied from decimeters to meters. Unfortunately, the lack of a validation dataset limited our ability to create a quantitative reliability map for this case. Figure 19 presents a qualitative sample map for illustrative purposes. The color scheme (light blue, yellow, red) on this map visually distinguishes between the classes (water, vegetation, bare soil), with each class corresponding to a specific data source. The accuracy in each area is directly linked to the reliability of the data source used.

6.4. Contributions and Limitations

The proposed approach provides several scientific advantages over the current state of the art, addressing key limitations in existing methodologies for riverbed modeling and morphodynamic analysis. A novel method capable of extrapolating a complete riverbed DTM was developed and successfully tested, offering a more detailed and accurate representation of the river’s morphology. This result represents a significant advancement in the field, enabling continuous monitoring and extraction of crucial morphodynamic data through a relatively straightforward and efficient image analysis process. This efficiency reduces the need for costly and time-consuming field measurements, making it accessible for a wider range of applications, including routine monitoring and large-scale river systems.

Figure 19. Qualitative example of a DTM reliability map, showing the different classes (water, bare soil, vegetation) with different color patterns, corresponding to specific data sources (RDB + water surface, photogrammetry, interpolation) characterized by different reliability. The DTM resulting from the proposed procedure is in the background.

Figure 19. Qualitative example of a DTM reliability map, showing the different classes (water, bare soil, vegetation) with different color patterns, corresponding to specific data sources (RDB + water surface, photogrammetry, interpolation) characterized by different reliability. The DTM resulting from the proposed procedure is in the background.

Moreover, the model proposed by Ragno et al. [30] was rigorously tested with morphodynamic data derived from the image analysis.

The overall agreement between the model and field observations is encouraging, as the model successfully captured the essential features of the bedform pattern, i.e., an alternating sequence of pools and riffles. The fact that the model predictions and field data exhibited comparable bar wavelengths in the initial kilometers suggests that the model, although simple, was capable of representing the main morphodynamic processes. The observed decrease in bar wavelength further upstream (beyond 2 km) was also captured by the model, although there were some discrepancies in the exact wavelength values.

A possible explanation for the larger differences in bar wavelengths between the model and the observations in the landward reach of the river could be related to the fact that, as shown in Figure 15, the river was not entirely straight and included several bends. While the seaward half of the river (from meter 0 to approximately meter 3600) was nearly straight, the bends were concentrated in the landward half (from meter 3600 to 4700). These bends significantly affected the morphodynamics of the river, potentially leading to the formation of different types of bedforms (the so called “forced bars”) that the model could not capture. These bars, unlike “free bars” that form spontaneously due to bed instability, are generated by the curvature of the river itself and typically do not migrate. Forced bars tend to remain stationary, influenced by the channel geometry, while free bars are mobile and can shift downstream over time. Additionally, free migrating bars are generally shorter than the meander wavelengths, which could explain the observed discrepancies in the model’s predicted bar wavelengths compared to real-world measurements in the bent section of the river. A key question is whether the bedforms in the landward half of the river reach are free bars, forced bars, or a combination of both. This distinction is important for understanding the dynamics of bedform evolution in the presence of bends. The fundamental research of [66,67] demonstrated that free and forced bars could coexist in rivers, provided that the meander sinuosity was not too pronounced. In channels with mild meanders, alternate free bars can migrate despite the bends, leading to complex interactions between these two types of bedforms. To clarify whether the bedforms in the landward half of the river reach are free, forced, or a combination of both, additional information is needed—particularly regarding their migration rates. Future research could resolve this issue by comparing successive satellite images of the river, which would allow us to track the movement (or lack thereof) of the bedforms over time. Such comparisons would provide a clearer understanding of the interaction between free and forced bars in the bent section of the river, offering valuable insights into the morphodynamic processes governing the river’s evolution.

Among the other simplifications included in the simulations, it is also worth mentioning the constant-channel-width assumption, which could also explain the slight deviations in bar wavelengths between the predicted and observed results. This hypothesis can be easily tested in the future by removing the constant-width constraint, as the model is already designed to accommodate channels with varying widths. Nevertheless, the alignment between the predicted and observed bar wavelengths reinforces the model’s utility in simulating riverbed morphodynamics.

The testing results indicate the model’s strong potential for accurately predicting the future response of the Roya River under various scenarios, including climate change-induced sea level rise and extreme discharge events. This is a critical breakthrough as it offers valuable insights into how river systems will adapt or respond to environmental changes, enabling better risk assessment and management strategies. Predicting future river responses is essential for sustainable water resource management, flood risk mitigation, and the preservation of riverine ecosystems, especially in the face of increasing climate variability. Additionally, the integration of this method with climate models opens the door for simulating a wider range of potential future scenarios, making it a versatile tool for both scientific research and practical decision-making.

The primary limitations of this study arise from the theoretical framework employed, which, due to its inherently linear nature, is unable to predict the equilibrium bar amplitude. The linear assumption, while effective for predicting early-stage instabilities in the flow and sediment transport, falls short in describing the full range of interactions that occur as the system evolves. Specifically, as the perturbation amplitude increases, nonlinear interactions between the flow field and the evolving bed topography become more prominent, complicating the behavior of the system. These nonlinear effects are critical in controlling the later stages of bar formation and growth, including the transition towards equilibrium states.

To overcome these limitations, the approach could benefit from incorporating nonlinear models that account for the complex interactions between flow dynamics and sediment transport. Nonlinear models, such as those based on perturbation theory or full numerical simulations, are capable of capturing the feedback mechanisms that emerge as bedforms grow and interact with the flow. These models would allow for the prediction of bar amplitudes and offer insights into the transient dynamics as the system evolves from an initial perturbed state towards equilibrium.

Further expansion of the study could investigate secondary processes such as sediment cohesion, the influence of vegetation, or variable flow regimes, which can significantly impact the equilibrium bar geometry. These factors, commonly found in natural river systems, can alter both the timescale and final state of morphodynamic adjustments. Extending the current framework to include these complexities would enhance the model’s predictive power, making it a more realistic tool for river management and restoration efforts.

Finally, one may wonder why using a simplified stability analysis instead of complete 3D commercial numerical code. In order to answer, it must be thus recalled that the present study represents the first phase of a larger project that aims to integrate a morphodynamic model with photogrammetry for improved river management. The photogrammetry tool offers continuous monitoring of the river, while the morphodynamic model supports management by providing quick and user-friendly simulations to address questions about bedform formation (e.g., if dredging occurs in one area, will deposition follow?), bedform migration (e.g., will the excavation shift over time? How quickly?), and long-term future scenarios (e.g., what changes can be expected if the river flow increases?). All these information are easily provided by simplified models. Once potential “fragility” areas are identified, the user would then need to turn to more detailed ad hoc simulations using three-dimensional commercial models to assess local morphodynamics. However, these models cannot provide the necessary preliminary insights within reasonable time frames or computational costs.

7. Conclusions

This study developed and validated a procedure for generating a seamless riverbed DTM by integrating various data sources, including an aerial multispectral orthophoto and a UAV photogrammetric elevation survey. This integration allowed for a comprehensive representation of the riverbed, encompassing both emerged and submerged sections, overcoming the limitations of individual data acquisition methods. The derived DTM proved highly useful, as evidenced by its application in extracting input parameters for a morphodynamic model and validating the model’s predictions against observed bar patterns. The research results indicated a qualitative alignment between the theoretical model and remote sensing data, underscoring the effectiveness of the integrated approach. This alignment suggests that combining these methods significantly improves our understanding of river morphology, hydrodynamics, and sediment dynamics. Such insights are crucial for informed decision-making in sustainable river management. The integrated approach not only enhances the accuracy of riverbed representations but also supports better management strategies by providing a comprehensive view of river dynamics and potential future changes.

Future research will focus on analyzing additional case studies and developing advanced remote sensing data analysis tools and methodologies to further enhance the understanding of river morphodynamics. The procedure for generating riverbed DTMs will be extended to other remote sensing imagery to assess river stability and bar migration. An important future objective is to validate the predictions of bar migration speed, here quantified as some cm/s using Ragno et al. [30]’s theoretical model, although this will depend on the availability of suitable data. The model will also be tested under different river discharge conditions and sea level elevations to evaluate bedform stability and river responses to various future scenarios. Furthermore, potential improvements to the model include incorporating vegetation growth over bars and accounting for varying bed roughness into the theoretical framework.