A Method for Riverway Erosion and Deposition Calculation Using Measured Cross-Sections

Abstract

The calculation of erosion and deposition in riverways plays a pivotal role in river morphology studies, comprehensive river management and flood safety. Some existing methods have certain limitations in terms of accuracy and applicability. To address these challenges, we propose a method for calculating riverway erosion and deposition based on measured cross-sectional terrain data. This method not only enables the calculation of changes in erosion and deposition along the riverway but also provides information on the spatial distribution of these changes. The validity of the proposed calculation method was assessed using measured bathymetric data. The results indicate a relative error of only 5.6% between the calculated and measured values for the total volume of erosion and deposition. A comparison with the results obtained using a cross-section method reveals that, with an average distance between adjacent sections of 1.0 km, the proposed method generally outperforms the cross-section method. The relative error in the total volume of erosion and deposition decreases from 19.2% with the cross-sectional method to 5.6% with our proposed method. When facing the need to calculate changes in riverway erosion and deposition, our approach offers a more accurate and flexible computational method.

1. Introduction

The calculation of erosion and deposition in riverways plays a pivotal role in studying river morphology, performing comprehensive river management and ensuring flood safety [1,2]. The erosion and deposition conditions are critical indicators for assessing the stability and safety of rivers, serving as essential bases for formulating river management and planning strategies [3,4,5].

Currently, there are numerous methods for calculating riverway erosion and deposition. For instance, Azarang et al. [6] simulated the erosion and deposition downstream of the Karkheh reservoir dam before construction using the MIKE11 model. Adib et al. [7] employed the Fluvial-12 model to compute erosion and deposition changes in the Karun River in Iran. However, the two most widely used methods are the sediment transport rate method and the cross-section method [8,9,10]. The sediment transport rate method involves calculating erosion and deposition based on measured sediment transport data at the inlet and outlet hydrological stations of a reach [11]. Due to the continuous nature of sediment data collection at hydrological stations, it can reflect the variations in erosion and deposition processes along the entire studied reach. However, the sediment transport rate method can only reveal the total volume of erosion and deposition over a long river segment between two hydrological stations [12]. It cannot capture along-riverway variations in erosion and deposition between the two stations. Moreover, the distance between two hydrological stations is typically long in actual rivers [13], making this method unsuitable for short river segments. Additionally, due to the current limitations in sediment transport data measurement technology, the measurement of sediment concentration is often influenced by the measurement point selection and sampling depth, resulting in low measurement accuracy and introducing errors into the calculation results [14]. Furthermore, these errors gradually accumulate over time. Importantly, this method, as it cannot account for changes in erosion and deposition caused by bedload transport, is only applicable to rivers where suspended sediment transport is predominant.

Simultaneously, in order to investigate the evolution of river erosion and deposition, it is common to set up fixed-observation cross-sections at regular distances along the river. With the terrain data from these observed cross-sections as a foundation, the cross-section method becomes a viable approach to calculating the erosion and deposition in a riverway. This method approximates the river reach between two measured cross-sections as a truncated cone. The areas of flow for the two cross-sections are treated as the two cross-sectional areas of the truncated cone. The distance between the two measured cross-sections serves as the length of the truncated cone. The volume of the truncated cone represents the channel storage volume [15,16]. The difference in the calculated channel storage volumes between the two instances corresponds to the erosion and deposition volume. For instance, Liu et al. [17] used the cross-section method to calculate the changes in riverbed erosion and deposition in the Lower Yellow River from 1965 to 2014 to study the impact of riverbed and floodplain changes on water and sediment conditions. Xia et al. [18] applied the cross-section method to calculate sediment erosion and deposition in the middle and lower reaches of the Yangtze River since the completion of the Three Gorges Project, supporting channel management and basin sediment management. Yuan et al. [19] used the cross-section method to study the trends in riverbed topographic changes in the Ningxia Plain Reaches of the Yellow River, calculating sediment erosion and deposition at different time periods and their variations. Lyu et al. [20] used the cross-section method to investigate riverbed changes in the middle reaches of the Yangtze River from Yichang to Chenglingji, providing insights into the impact of the Three Gorges Project on downstream riverbed morphodynamic adjustments. However, through the analysis of the aforementioned research outcomes, we identified some shortcomings in employing the cross-section method for riverway erosion and deposition calculations. Firstly, it fails to reflect the meandering morphology of rivers since it directly connects two cross-sections with a straight line when conceptualizing the truncated cone. Natural rivers are typically not straight, leading to potential calculation errors [21,22]. Secondly, it does not depict the spatial distribution of erosion and deposition changes along the calculated riverway. The cross-section method calculates the erosion and deposition between adjacent cross-sections and accumulates these values to obtain the total erosion and deposition volume for the entire riverway [23]. However, understanding the distribution of erosion and deposition is crucial for studying river channel changes. This can provide a basis for targeted management measures for river managers.

Based on this, this paper proposes a method for calculating riverway erosion and deposition based on measured cross-sectional terrain data. This method is designed for the computation of riverway erosion and deposition. It not only considers the meandering morphology of rivers by introducing auxiliary lines, such as river boundaries and water edges, but also incorporates a quadrilateral grid to obtain the spatial distribution of erosion and deposition changes in the riverway. In comparison to the cross-section method, this approach can provide more accurate calculation results and offer information about the spatial distribution of erosion and deposition changes that the cross-section method cannot provide.

2. Riverway Erosion and Deposition Calculation Method

For a river segment requiring the calculation of erosion and deposition changes, based on cross-sectional terrain data measured at two different times, the following steps can be employed to calculate the erosion and deposition variations (designated as MT₁ and MT₂ for the first and second measurement times, respectively). Firstly, the riverway boundaries and water edges corresponding to MT₁ and MT₂ are utilized as auxiliary lines to create two sets of quadrilateral grids. Subsequently, the elevation information for all nodes on these two grids is calculated using the cross-sectional terrain data obtained at MT₁ and MT₂. Finally, the erosion and deposition conditions of the target river segment during the MT₁–MT₂ period are determined using the two sets of grids corresponding to MT₁ and MT₂. The elevation of any point within the riverbed is influenced by the measured cross-sectional terrain data of its two adjacent measured sections. The entire river segment can be subdivided into multiple sub-segments using the observed cross-sections. These sub-segments share a consistent calculation method, illustrating that one sub-segment serves as an example. The detailed process can be roughly divided into the following five stages (see Figure 1):

Stage 1—Data preparation: The planar coordinate data for the boundaries on both sides of the riverway corresponding to MT₁ and MT₂, as well as the planar coordinate data for the left and right water edges corresponding to MT₁ and MT₂ and the measured cross-sectional terrain data corresponding to MT₁ and MT₂, must be prepared in advance.

Stage 2—Thalweg generation: The planar coordinate data for the left and right water edges corresponding to MT₁ and MT₂ and the measured cross-sectional terrain data corresponding to MT₁ and MT₂ are used to calculate the thalweg for the respective measurement times.

Stage 3—Riverway grid division: By incorporating the riverway boundaries, water edges and thalweg as auxiliary lines, quadrilateral grids corresponding to MT₁ are generated on the target river segment, effectively reflecting the course of the riverway. Similarly, quadrilateral grids corresponding to MT₂ are obtained.

Stage 4—Riverway grid node elevation interpolation: Measured cross-sectional terrain data corresponding to MT₁ and MT₂ are interpolated onto the grid nodes located on the measured cross-sections using distance-weighted interpolation methods. Subsequently, distance-weighted interpolation methods are employed in the water flow direction to obtain elevation information for all grid nodes between the two measured cross-sections.

Stage 5—Riverway erosion and deposition calculation: The planar coordinates of the two sets of grids are matched and calculated with node elevation data corresponding to MT₁ and MT₂. The erosion and deposition volume and spatial distribution between the two measured cross-sections during the MT₁–MT₂ period are determined.

2.1. Data Preparation

Before commencing the calculations, the cross-sectional terrain data corresponding to MT₁ and MT₂, the planar coordinate data for the boundaries on both sides of the riverway corresponding to MT₁ and MT₂, and the planar coordinate data for the water edges on both sides of the riverway corresponding to MT₁ and MT₂ need to be prepared.

2.2. Thalweg Generation

In order to make the generated quadrilateral grids better adapt to the course of the river and changes in the terrain, water edges and the thalweg are introduced as auxiliary lines in the grid generation process, in addition to the riverway boundaries. The riverway boundaries and water edges can be obtained through field measurements or by interpreting remote sensing images [24,25,26,27]. In contrast, obtaining the thalweg can be more challenging. Thalwegs are submerged beneath the water and are generally difficult to directly obtain. Field measurements of thalwegs can be labor-intensive and resource-intensive. Here, we propose a method for determining the thalweg using water edges and measured cross-sectional terrain data, with an example focused on the subsection between measured cross-sections c and c + 1. Furthermore, in accordance with the conventions of river hydrodynamics, we define the left side as the left bank and the right side as the right bank when facing the water flow; this convention will be consistently followed.

The left and right water edges are equally divided into N segments by length [28], resulting in nodes along the left and right water edges, denoted by $(x_{LW}(i),y_{LW}(i))$ and $(x_{RW}(i),y_{RW}(i))$, where i = 1,2,3,…, N + 1 (see Figure 2a). In fact, the value of N should be appropriately selected based on the distance between the two measured cross-sections. A larger N results in a smoother thalweg, but it also leads to a higher computational workload. Points at the corresponding positions of the left and right water edges are sequentially connected to form N + 1 subsections (see Figure 2b). Subsequently, a dimensionless parameter, denoted by T, is introduced to represent the deviation of the thalweg point between the two water edges, calculated using Equation (1). The T values range between 0 and 1, with values closer to 0 indicating that the thalweg point is closer to the left water edge and values closer to 1 indicating proximity to the right water edge. The T values, T(c) and T(c + 1), are computed separately for cross-sections c and c + 1.

$$T \hspace{0.17em}=\hspace{0.17em}(L_T-L_{LW})/(L_{RW}-L_{LW})$$

(1)

where $L_T$ represents the horizontal distance of the thalweg point from the left endpoint of the observed cross-section, $L_{LW}$ denotes the horizontal distance of the left water edge from the left endpoint of the observed cross-section, and $L_{RW}$ represents the horizontal distance of the right water edge from the left endpoint of the observed cross-section.

Subsequently, interpolation is employed according to Equation (2) to calculate the T values, T(i), for all N + 1 subsections between the two observed cross-sections. Based on the T values for each subsection and the planar coordinates, $(x_{LW}(i),y_{LW}(i))$ and $(x_{RW}(i),y_{RW}(i))$, of their left and right endpoints, the planar coordinates of each thalweg point, denoted by $(x_T(i),y_T(i))$, are computed using Equations (3) and (4), where i = 1,2,3,…, N + 1. Subsequently, the generated thalweg points are sequentially connected to form the thalweg, as illustrated in Figure 2c.

$$T(i) \hspace{0.17em}=\hspace{0.17em}T(c)+\frac{i-1}{N}\cdot (T(c+1)-T(c))$$

(2)

$$x_T(i)=x_{LW}(i)+T(i)\cdot (x_{RW}(i)-x_{LW}(i))$$

(3)

$$y_T(i)=y_{LW}(i)+T(i)\cdot (y_{RW}(i)-y_{LW}(i))$$

(4)

2.3. Riverway Grid Division

In the process of subdividing the target riverway using quadrilateral grids, auxiliary lines such as the riverway boundaries, water edges and thalweg are introduced, assuming that there are a total of M such auxiliary lines. Each of these auxiliary lines is equally divided into N segments by length. The coordinates of the ith node on the mth auxiliary line are denoted by $(x_A(i,m),y_A(i,m))$, where i = 1, 2, 3, …, N + 1, and m = 1, 2, 3, …, M (see Figure 3a). Subsequently, nodes with the same index on adjacent auxiliary lines, $(x_A(i,m),y_A(i,m))$ and $(x_A(i,m+1),y_A(i,m+1))$, are connected in sequence, resulting in a preliminary grid, as shown in Figure 3b. It is evident that a grid subdivided in this manner exhibits relatively large spacing between the two auxiliary lines, which significantly impacts the accuracy of the computation results. Therefore, refinement is required in the region between the two auxiliary lines. Here, we divide the lines connecting corresponding nodes on every two auxiliary lines into Q segments. Of course, the number of segments between different pairs of auxiliary lines may vary and can be adjusted based on actual distances to ensure a uniform grid density across the entire region. The new grid nodes $(x_{NA}(i,q,m)\hspace{0.17em},y_{NA}(i,q,m)\hspace{0.17em})$ between two auxiliary lines can be calculated using Equations (5) and (6), where i = 1, 2, 3, …, N + 1; q = 1, 2, 3, …, Q − 1; and m = 1, 2, 3, …, M − 1 (see Figure 3c). For ease of programming, the arrays $x_{NA}(i,q,m)\hspace{0.17em}$ and $y_{NA}(i,q,m)\hspace{0.17em}$ are reduced in dimensionality. The nodes on the auxiliary lines $(x_A(i,m),y_A(i,m))$ are also included. This results in a restructured set of coordinates for all grid nodes, denoted by $(x_{node}(i,j),\hspace{0.17em}{y}_{node}(i,j))$, as shown in Figure 3c, where i = 1, 2, 3, …, N + 1; j = 1, 2, 3, …, J; and J = 1+${\displaystyle \sum Q(m)}$.

$$x_{NA}(i,q,m)\hspace{0.17em}=\hspace{0.17em}{x}_A(i,m)+(q/Q)\times (x_A(i,m+1)-x_A(i,m))$$

(5)

$$y_{NA}(i,q,m)\hspace{0.17em}=\hspace{0.17em}{y}_A(i,m)+(q/Q)\times (y_A(i,m+1)-y_A(i,m))$$

(6)

2.4. Riverway Grid Node Elevation Interpolation

Interpolating the actual cross-sectional terrain data to the grid nodes involves two steps: firstly, interpolating the actual cross-sectional terrain data to the grid nodes located on the observed cross-sections; secondly, using the distance between adjacent observed cross-sections as weights and calculating the elevations for all grid nodes between two observed cross-sections. The specific steps are as follows:

Step 1: Elevation interpolation for grid nodes located on the observed cross-sections. For each grid node $(x_{node}(0,j),y_{node}(0,j))$ located on the observed cross-section c, its horizontal distance $d_c(j)$ to the left endpoint $\hspace{0.17em}(xL(c),yL(c))$ can be calculated using Equation (7). Similarly, for each grid node $(x_{node}(N+1,j),y_{node}(N+1,j))$ located on the observed cross-section c + 1, its horizontal distance $d_{c+1}(j)$ to the left endpoint $\hspace{0.17em}(xL(c+1),yL(c+1))$ can be calculated using Equation (8). For any grid nodes located on cross-section c, if $d_c(j)$ falls between the horizontal distances of the two observed points on cross-section c, namely, $L_c(r)<=d_c(j)<L_c(r+1)$, the elevation $z_c(j)$ can be calculated using Equation (9) according to the two corresponding measured elevation values $z_m(r)$ and $z_m(r+1)$ (see Figure 4). Similarly, the elevation $z_{c+1}(j)$ can be calculated for grid nodes located on cross-section c + 1.

$$d_c(j)\hspace{0.17em}=\hspace{0.17em}\sqrt{{(x_{node}(0,j)-xL(c))}^2+{(y_{node}(0,j)-yL(c))}^2}\hspace{0.17em}$$

(7)

$$d_{c+1}(j)\hspace{0.17em}=\hspace{0.17em}\sqrt{{(x_{node}(N+1,j)-xL(c+1))}^2+{(y_{node}(N+1,j)-yL(c+1))}^2}\hspace{0.17em}$$

(8)

$$z_c(j)\hspace{0.17em}=\hspace{0.17em}{z}_m(r)+\frac{d_c(j)-L(r)}{L(r+1)-L(r)}\hspace{0.17em}\times \left(z_m(r+1)-z_m(r)\right)$$

(9)

Step 2: Elevation interpolation for grid nodes between adjacent observed cross-sections. The cumulative distances from each grid node located on every longitudinal (in the flow direction) grid line to cross-section c are computed. The cumulative distance for the first grid node on each longitudinal grid line to cross-section c is set to zero, while for the subsequent nodes, it is calculated using Equation (10). In essence, the cumulative distance $d(N+1,j)$ for the last grid node on each longitudinal grid line is equivalent to the total length of that particular longitudinal grid line (Figure 5). Subsequently, the cumulative distances of the target interpolation grid nodes to both observed cross-sections are employed as weights, using Equation (11) to interpolate and calculate their elevation values $z(i,j)$.

$$d(i,j)\hspace{0.17em}=\hspace{0.17em}d(i-1,j)+\sqrt{{(x_{node}(i,j)-x_{node}(i-1,j))}^2+{(y_{node}(i,j)-y_{node}(i-1,j))}^2}\hspace{0.17em}$$

(10)

$$z(i,j)\hspace{0.17em}=\hspace{0.17em}{z}_c(j)+\frac{(d(i,j)}{d(N+1,j))}\cdot (z_{c+1}(j)-z_c(j))$$

(11)

2.5. Riverway Erosion and Deposition Calculation

Using the aforementioned methods, we obtain the coordinates of the grid nodes, with $(x^{\prime }(i,j),y^{\prime }(i,j),z^{\prime }(i,j))$ corresponding to MT₁ and $(x(i,j),y(i,j),z(i,j))$ corresponding to MT₂. Due to the differing positions of the water edges for MT₁ and MT₂, the location of the thalweg naturally exhibits some variation. While the generated grids have the same number of rows and columns, the planar coordinate positions of the grid nodes differ somewhat. To calculate the elevation change for grid nodes corresponding to MT₁ and MT₂, it is necessary to align the planar coordinates of the nodes from these two sets of grids by means of interpolation. This enables the direct calculation of the elevation difference for the corresponding nodes. The specific procedure involves interpolating the grid nodes $\left(x^{\prime }(i,j),y^{\prime }(i,j),z^{\prime }(i,j)\right)$ corresponding to MT₁ onto the grid nodes $\left(x(i,j),y(i,j)\right)$ corresponding to MT₂, yielding the grid node coordinates $\left(x(i,j),y(i,j),z^{″}(i,j)\right)$. A four-point natural neighbor interpolation method is employed for this purpose. In essence, this method selects the four closest points, $\left(x_1,y_1,z_1\right)$, $\left(x_2,y_2,z_2\right)$, $\left(x_3,y_3,z_3\right)$ and $\left(x_4,y_4,z_4\right)$, from the grid nodes $\left(x^{\prime }(i,j),y^{\prime }(i,j),z^{\prime }(i,j)\right)$ corresponding to MT₁ to the target point $\left(x(i,j),y(i,j)\right)$ based on their planar Euclidean distances $d_1$, $d_2$, $d_3$ and $d_4$ to the target point. Subsequently, the elevation $z″(i,j)$ at the target point can be calculated as per Equation (12). Furthermore, the change in riverbed elevation $D_z(i,j)$ at each grid node can also be calculated using Equation (13). If $D_z(i,j)<0$, the node exhibits erosion; if $D_z(i,j)>0$, the node indicates deposition; if $D_z(i,j)=0$, the node remains unchanged. The distribution of erosion and deposition in the entire sub-river segment is illustrated in Figure 6.

$$z″(i,j)\hspace{0.17em}={\displaystyle \sum _{k=1}^4\frac{z_k}{d_k}}/{\displaystyle \sum _{k=1}^4\frac{1}{d_k}}$$

(12)

$$D_z(i,j)\hspace{0.17em}=\hspace{0.17em}z(i,j)-z″(i,j)$$

(13)

As the plane coordinates of the grid nodes corresponding to MT₁ and MT₂ have already been aligned, the calculation element for sediment erosion and deposition, denoted by ABCDEFGH in Figure 7, is defined by its eight vertices, whose coordinates are listed in

Table 1. Coordinates of the vertices on the calculation element.

Vertices	X	Y	Z
A	$x(i,j)$	$y(i,j)$	$z(i,j)$
B	$x(i,j+1)$	$y(i,j+1)$	$z(i,j+1)$
C	$x(i+1,j+1)$	$y(i+1,j+1)$	$z(i+1,j+1)$
D	$x(i+1,j)$	$y(i+1,j)$	$z(i+1,j)$
E	$x(i,j)$	$y(i,j)$	$z″(i,j)$
F	$x(i,j+1)$	$y(i,j+1)$	$z″(i,j+1)$
G	$x(i+1,j+1)$	$y(i+1,j+1)$	$z″(i+1,j+1)$
H	$x(i+1,j)$	$y(i+1,j)$	$z″(i+1,j)$

. This unit represents a polyhedron with four vertices, A, B, C and D, corresponding to the bed elevations of MT₂, and four vertices, E, F, G and H, representing the bed elevations of MT₁. The volume of this polyhedron, denoted by ABCDEFGH, represents the sediment volume change within this calculation element. Since the vertices A, B, C and D may not lie in the same plane, the vertices E, F, G and H may not be coplanar, and it is necessary to further divide this polyhedron into two triangular prisms, ABCEFG and ACDEGH. Taking the polyhedron ABCEFG as an example, the upper-surface vertices A, B and C are reordered according to their node elevations in ascending order, resulting in the sequence $a\hspace{0.17em}(x_a,y_a,z_a)$, $b\hspace{0.17em}(x_b,y_b,z_b)$ and $c\hspace{0.17em}(x_c,y_c,z_c)$, ensuring that $z_a<z_b<z_c$. Likewise, the lower-surface vertices E, F and G are reordered based on their elevations in ascending order, resulting in the sequence $e\hspace{0.17em}(x_e,y_e,z_e)$, $f\hspace{0.17em}(x_f,y_f,z_f)$ and $g\hspace{0.17em}(x_g,y_g,z_g)$, ensuring that $z_e<z_f<z_g$. The volume $V_1$ of this polyhedron can be calculated according to Equations (14) and (17). Similarly, the volume $V_2$ of the polyhedron ACDEGH is determined. The erosion and deposition volume $V$ within the sub-river segment between the two observation cross-sections is computed as per Equation (18).

$$S=(1/2)\cdot \left|x_a\cdot (y_b-y_c)+y_b\cdot (y_c-y_a)+y_c\cdot (y_a-y_b)\right|$$

(14)

$$V_{abc}=(z_a\cdot S)+(1/6)\cdot (z_b+z_c-2z_a)\cdot \left|(y_b-y_c)\cdot x_a-(x_b-x_c)\cdot y_a+x_b\cdot y_c-y_b\cdot x_c\right|$$

(15)

$$V_{efg}=(z_e\cdot S)+(1/6)\cdot (z_f+z_g-2z_e)\cdot \left|(y_f-y_g)\cdot x_e-(x_f-x_g)\cdot y_e+x_f\cdot y_g-y_f\cdot x_g\right|$$

(16)

$$V_1=V_{abc}-V_{efg}$$

(17)

$$V={\displaystyle \sum _{i=1,j=1}^{N,J-1}\left(V_1(i,j)+V_2(i,j)\right)}$$

(18)

where S is the projected area of the upper surfaces ABC and lower surfaces EFG onto the horizontal reference plane.

3. Results

To validate the rationality and computational accuracy of the proposed method, we selected a 13.8 km reach of the Pearl River in China for verification. The reason for choosing this particular reach is that we have measured bathymetric data for this reach at two different measurement times. Figure 8a,b, respectively, illustrate the measured bathymetric data for this reach on 26 January 2012 and 15 February 2020. Using an average distance between adjacent sections (ADAS) of approximately 1.0 km, we extracted 15 cross-sections from each of these two datasets. These can be considered the measured cross-sectional terrain data corresponding to the two measurement times. These two sets of cross-sections will be used to calculate the erosion and deposition changes in this reach from 26 January 2012 to 15 February 2020. Figure 8a,b also include the riverway boundaries and water edges corresponding to the measurement times of the bathymetric data. Based on the extracted cross-sectional terrain data, the results obtained using the proposed method will be compared with the measured erosion and deposition changes, thereby verifying the rationality and computational accuracy of the method.

Figure 9 presents a comparison between the calculated values and the measured values of erosion and deposition changes in the riverway during the specified period. In Figure 9, it can be observed that the calculated values generally align well with the measured values for overall erosion and deposition changes in the reach. The calculated total volume of erosion and deposition for this reach is 2.625 × 10⁶ m³, while the measured value is 2.485 × 10⁶ m³, resulting in a relative error of only 5.6%. For specific segments of the reach, between cross-sections 1 to 5 and 10 to 15 (in relatively straight sections), the agreement between calculated and measured values is higher. The regions of erosion and deposition changes are consistent. However, between cross-sections 4 and 5, there is a significant area with erosion and deposition changes in the measured values (Figure 9b) that does not appear in the distribution map obtained using the proposed method (Figure 9a). This discrepancy is attributed to the fact that, in the method proposed in this paper, the riverbed terrain is interpolated based on the terrain data of adjacent cross-sections. Our interpolated terrain is always controlled by the elevation of the measured cross-sections. Natural riverway changes in bed erosion and deposition are influenced by both water and sediment conditions and bed boundary conditions. Cross-sections cannot sample every specific shape of the riverbed, leading to the neglect of some terrain changes between cross-sections. This can also be observed in Figure 9c, where the differences between the measured values and calculated values primarily occur in the region between two adjacent cross-sections. In the vicinity of the cross-sections, such differences are minimal. In the segment between cross-sections 5 and 10 (with a river center continent), the calculated and measured volume of the erosion and deposition are 0.016 × 10⁶ m³ and 0.035 × 10⁶ m³, respectively. However, due to the presence of a river center continent in this reach, the meandering morphology of the riverway is more complex. The spatial distribution of erosion and deposition changes obtained by the interpolation calculation differs somewhat from that calculated from measured values. Figure 10 provides a comparison between the calculated and measured volumes of erosion and deposition between adjacent cross-sections. To evaluate the average level of proximity between the calculated values and measured values for the erosion and deposition volume in various river segments, the Nash efficiency coefficient (NSE) was calculated, yielding a value of 0.88. This indicates that the proposed method in this paper exhibits good computational quality. Overall, the trend of changes is generally consistent, but there are occasional significant differences between certain pairs of adjacent cross-sections, such as between cross-sections 4 and 5, where interpolation fails to capture a significant area of erosion and deposition changes, resulting in a calculated value that is greater than the measured value.

4. Discussion

4.1. Comparison with the Cross-Section Method

We also employed the cross-section method proposed by Yang et al. [16] to calculate the erosion and deposition changes in the reach mentioned in Section 3, obtaining the total erosion and deposition volume for the entire reach (

Table 2. Total erosion and deposition volume using the proposed method and the cross-section method.

Method	Total Erosion and Deposition Volume
	Calculated Value (106 m3)	Measured Value (106 m3)	Absolute Error 1 (106 m3)	Relative Error 2 (%)
Cross-section method	2.961	2.485	0.476	19.2
Method in this paper	2.625		0.140	5.6

¹ Absolute Error = Calculated Value − Measured Value. ² Relative Error = (Calculated Value − Measured Value)/Measured Value.

). In terms of the total erosion and deposition volume in the reach, the cross-section method calculated a value of 2.961 × 10⁶ m³, while the proposed method in this paper calculated a value of 2.625 × 10⁶ m³. These values represent a 19.2% and 5.6% overestimation, respectively, compared to the measured value of 2.485 × 10⁶ m³. This indicates that the proposed method in this paper exhibits a significant improvement in computational accuracy compared to the cross-section method. Figure 11 provides a comparison of erosion and deposition volume between adjacent cross-sections using different calculation methods. In Figure 11, it can be observed that the NSE for the proposed method in this paper is 0.88, while the NSE for the cross-section method is 0.71. This suggests that the proposed calculation method in this paper is superior to the cross-section method, with a higher level of agreement between the calculated values of erosion and deposition volume between adjacent cross-sections and the measured values. In relatively straight river sections, both the proposed method in this paper and the cross-section method yield similar results that align well with the measured values. However, in river sections with significant changes in the meandering morphology, the results obtained by the proposed method in this paper are closer to the measured values compared to the results from the cross-section method. This is because the cross-section method calculates erosion and deposition changes by approximating the river segment between two cross-sections as a truncated cone, relying heavily on the representative nature of the cross-sections and neglecting variations in width along the reach and the curvature of the river channel. In contrast, the method proposed in this paper divides the entire reach into grids, with the grid orientation controlled by auxiliary lines such as the riverway boundaries, fully considering changes in the channel width and channel curvature along the reach. In this respect, it is superior to the cross-section method.

4.2. Grid Size Sensitivity Analysis

In a grid-based numerical method, the grid size can sometimes affect the accuracy of the results. In order to understand the influence of grid size on the calculation results in the proposed method, seven different grid sizes were set, and their corresponding calculation results were compared to analyze the grid size sensitivity. In the case study in Section 3 of this paper, the measured average spacing of the cross-sectional terrain data is 10.0 m. Based on this, grids with average sizes of 80.0 m, 40.0 m, 20.0 m, 10.0 m, 5.0 m, 3.0 m and 1.0 m were generated, and their corresponding erosion and deposition volumes were calculated (

Table 3. Total erosion and deposition volume with different grid sizes.

The Average Grid Size (m)	Total Erosion and Deposition Volume
	Calculated Value (106 m3)	Measured Value (106 m3)	Absolute Error 1 (106 m3)	Relative Error 2 (%)
80	3.021	2.485	0.536	21.6
40	2.886		0.401	16.1
20	2.715		0.230	9.3
10	2.624		0.139	5.6
5	2.625		0.140	5.6
3	2.625		0.140	5.6
1	2.625		0.140	5.6

¹ Absolute Error = Calculated Value − Measured Value. ² Relative Error = (Calculated Value − Measured Value)/Measured Value.

). It can be observed that as the grid size decreases, the error between the calculated values and the measured values decreases. However, when the grid size is less than 10.0 m, a further reduction in grid size does not result in significant changes in the outcomes. This is because, in the calculation of erosion and deposition according to the proposed method in this paper, all results are derived from cross-sectional terrain data. When the grid size is larger than the distance between adjacent measurement points on the cross-section, there is a certain degree of generalization of the cross-sectional terrain data, leading to a decrease in data resolution and introducing errors into the calculation results. When the grid size is smaller than this distance, a limited number of data points on the cross-section are still used as source data, but better results are not achieved. Thus, when using this method, the optimal grid size should be consistent with the measured average spacing of cross-sectional terrain data.

4.3. Impact of the Distance between Adjacent Cross-Sections on Calculation Accuracy

The distance between adjacent cross-sections is one of the crucial factors influencing the calculation accuracy. To compare the effect of ADAS on the accuracy of the proposed calculation method in this paper, the ADAS used in Section 3 was successively reduced from 1.0 km to 0.8 km, 0.6 km and 0.4 km. Correspondingly, the number of cross-sections increased from 15 to 19, 24 and 36, respectively.

Table 4. Total erosion and deposition volume with different ADASs.

Number of Cross-Sections	ADAS (km)	Total Erosion and Deposition Volume
		Calculated Value (106 m3)	Measured Value (106 m3)	Absolute Error 1 (106 m3)	Relative Error 2 (%)
15	1.0	2.625	2.485	0.140	5.6
19	0.8	2.605		0.120	4.8
24	0.6	2.397		−0.088	−3.5
36	0.4	2.423		−0.062	−2.5

¹ Absolute Error = Calculated Value − Measured Value. ² Relative Error = (Calculated Value − Measured Value)/Measured Value.

provides the calculated erosion and deposition volume for the reach using different ADAS. It can be observed that as the ADAS decreases (i.e., the number of cross-sections increases), the calculated changes in erosion and deposition volume gradually approach the measured values. When the ADAS is 1.0 km, 0.8 km, 0.6 km and 0.4 km, the calculated changes in reach erosion and deposition are 2.625 × 10⁶ m³, 2.605 × 10⁶ m³, 2.397 × 10⁶ m³ and 2.423 × 10⁶ m³, respectively. These represent deviations of 5.6%, 4.8%, −3.5% and −2.5% from the measured value of 2.485 × 10⁶ m³. Figure 12 presents a comparison between the calculated and measured values of erosion and deposition volume between adjacent cross-sections for different ADASs. It is evident that as the ADAS decreases, the degree of agreement between the calculated values and the measured values increases. The ADAS decreases from 1.0 km to 0.8 km, 0.6 km and 0.4 km, resulting in an increase in the NES from 0.88 to 0.90, 0.92 and 0.95, respectively. Figure 13 illustrates the erosion and deposition distribution obtained using different ADASs. It can be observed that as the ADAS decreases, the degree of agreement between the calculated values and the measured values of erosion and deposition distribution continually improves. This allows for a more detailed representation of changes in erosion and deposition.

Therefore, it can be concluded that when calculating riverway erosion and deposition using the proposed method in this paper, the results are already superior to those obtained through the cross-section method when ADAS = 1.0 km. The results obtained with smaller ADASs are even closer to the measured values. In practical applications of this method, an appropriate ADAS should be chosen based on the characteristics of the river and the computational distance. With the development of river digital twinning and smart water management, more accurate data on river erosion and deposition changes are needed as support. When economic, human and time costs permit, using high-resolution river bathymetric data to reflect changes in riverway erosion and deposition is a preferable choice. However, if only cross-sectional terrain data are available, the method proposed in this paper remains a good approach to obtaining accurate erosion and deposition changes and their spatial distribution.

5. Conclusions

In this paper, a method for calculating riverway erosion and deposition using measured cross-sections is proposed, and detailed calculation steps are provided. The method generates a quadrilateral grid in the calculation domain by introducing auxiliary lines, such as riverway boundaries and water edges. It not only calculates accurate erosion and deposition changes but also provides their spatial distribution.

The proposed method was validated using measured data. The results show that the relative error between the calculated and measured values of the total erosion and deposition volume is only 5.6%. A comparison with the results of the cross-section method reveals that, at a cross-section spacing of 1.0 km, the overall accuracy of the proposed method is better than that of the cross-section method. The impact of the cross-section spacing on the results of the proposed method is discussed, showing that as the distance between adjacent cross-sections decreases from 1.0 km to 0.4 km, the relative error between the calculated and measured values of total erosion and deposition volume decreases from 5.6% to 2.5%.

When faced with the need to calculate changes in riverway erosion and deposition, the proposed method based on measured cross-sections can provide accurate data on erosion and deposition changes and their spatial distribution.