The initial mesh generation divides the computation domain into two regions. The first is the channel region, where the mesh has a higher spatial resolution since it is where key morphodynamic processes occur. The maximum size of the triangular elements in the channel region, d _s is defined as d _s = b/5, where b is the river channel half-width. The second region is the floodplain region, where a coarser resolution is used, with a size of triangular elements, d _f between 3d _s and 5d _s to reduce computational effort. The mesh size transition between the two regions is progressively increased to minimize numerical instability.

The initial river bathymetry and floodplain topography are mapped onto the initial mesh within the framework of a Cartesian coordinate system. In cases where digital elevation model (DEM) data for the floodplain is not accessible, an alternative method for defining the valley slope, S _v , is based on the equation S _v = S _c ⋅Ω, where S _c is the longitudinal channel slope and Ω is the channel sinuosity. The initial demarcation of the river bank lines is required to separate the channel and floodplain domains. These bank lines are subject to dynamic modifications enacted by the channel migration submodel.

TELEMAC-2D and SISYPHE are mature open-source models developed by EDF-R&D (France). More details about TELEMAC-2D and SISYPHE can be found in the two manuals by Lang (2013) and Tassi and Villaret (2014), respectively. We will briefly introduce the governing equations of the two models below.

TELEMAC-2D solves the two-dimensional time-dependent shallow water equations (Eqs. 1–3) using the h-type finite element method (Lang, 2013) for spatial discretization. The mass and momentum conservation equations are: ∂ h ∂ t + U ⋅ ∇ h + h ∇ ⋅ U = S r c h (1) ∂ u ∂ t + U ⋅ ∇ u = S r c x − g ∂ H ∂ x + 1 h ∇ ⋅ h v t ∇ u (2) ∂ v ∂ t + U ⋅ ∇ v = S r c y − g ∂ H ∂ y + 1 h ∇ ⋅ h v t ∇ v (3) where h is the water depth; t is the time; U = (u, v) is the depth-averaged velocity vector in Cartesian (x- and y-direction) coordinates, respectively; Src _h , Src _x and Src _y are the source/sink terms in mass and momentum conservation (bottom roughness is wrapped into the two terms Src _x and Src _y ); g is the gravitational acceleration; H is the water surface elevation; ν _t is the turbulent eddy viscosity. We used a k-ϵ model averaged over the vertical for turbulence closure in the 2D model utilized in this work (Hervouet, 2007).

SISYPHE computes the sediment transport (as bedload and suspended load for non-cohesive or cohesive sediments) and the evolution of the riverbed (Tassi and Villaret, 2014). Among several classical bedload transport formulae, the Meyer-Peter-M u ̈ ller formula (Eq. 4) was selected for validation of the cases analyzed herein. Φ b = 0 i f θ ′ ⩽ 0.047 8 θ ′ − 0.047 3 / 2 i f θ ′ > 0.047 (4) where Φ _b is the non-dimensional sand transport rate and θ′ is the Shields parameter.

The Exner equation (Eq. 5) is used in the SISYPHE model to solve bed evolution: ∂ z b ∂ t = − 1 1 − λ ∇ ⋅ Q b (5) where z _b is the bed elevation; t is the time; λ is the material porosity; Q _b = (Q _bx , Q _by ) is the vector of unit sediment discharges in x- and y-direction, respectively.

The formulation of Soulsby (1997) and Talmon et al. (1995) was used to correct the magnitude and the direction of the bed load transport, respectively, and Engelund’s formula to evaluate the deviation of sediment transport and the depth-averaged flow due to secondary flow.

The computation of the flow field is performed by TELEMAC-2D and the computation of Exner equation is addressed by SISYPHE. Flow field and bed morphology have different time scales; hydraulic perturbations occur faster than the bed evolution processes. For that reason, the computation of flow and bed evolution can be decoupled. Once the flow variables are solved, they are passed to SISYPHE to solve bedload and bed evolution.

Planform channel migration occurs by fluvial erosion and bank geotechnical processes causing bank retreat, and by depositional processes then produce bank advance (Parker et al., 2011). Here MEANDRE is used (Langendoen et al., 2016) to solve bank erosion and accretion as expressed by Eqs. 6, 7: E = M e τ τ c e − 1 i f τ > τ c e 0 i f τ ⩽ τ c e (6) A = M a 1 − τ τ c a i f τ < τ c a 0 i f τ ⩾ τ c a (7) where E and A are the bank erosion and accretion rates [m/s], respectively; M _e and M _a are the erosion-rate and accretion-rate coefficient (m/s), respectively; τ is the shear stress exerted by the flow on the banks; τ _ce is the critical shear stress for erosion of the material on the banks, on the other hand, τ _ca is a conceptual parameter equivalent to critical shear stress, on the contrary applied for accretion. The relation between the two previous critical shear stresses is important since their net effect is the widening or narrowing of the channel. Typically, we let M _e = M _a and τ _ce = τ _ca to maintain a dynamically uniform channel width. M _e , M _a , τ _ce , and τ _ca are model calibration parameters in MEANDRE, and they are typically chosen to match the observed channel migration rates in the study river through an analysis of multi-year remotely-sensed images.

Besides the flow field and bed evolution, bank migration (planform dynamics) is involved in the morphological evolution of rivers and it has the largest time scale. The algorithm proposed for the erosion and accretion processes also has a different coupling period from bed morphology and flow field computations. During the step of the computation of the erosion and accretion rates, we can obtain the migration distances of the river through Eq. 8. The migration of the banks determined by d, triggers adaptation of the elements of the mesh and reconfiguration of its topology (adaptation of nodes and reconfiguration of elements of the mesh). d = E ⋅ Δ t ⋅ e n e r o s i o n A ⋅ Δ t ⋅ e n a c c r e t i o n (8) where d is the migration vector; Δt is time step in the model; and e _n is the unit vector normal to the bank.

Previous studies have proposed an algorithm for detecting imminent cutoffs in numerical modeling (e.g., Sun et al., 1996; Camporeale et al., 2005). In this study, we utilized a modified version of Camporeale et al. (2005)’s algorithm, which uses the outer-bank line instead of the channel centerline and optimizes the sweeping range. These modifications significantly reduce iteration times compared to enumerating all nodes to calculate their distances and locate the minimum distance that identifies an imminent cutoff. The algorithm is briefly described below (Figure 2):

1. The outer-bank lines are discretized at equal length intervals (d _s in Figure 2) and n nodes that are sorted onto a uniform square grid board. The size of each cell is exactly equal to the threshold distance of cutoff d _c , where d _c is set to be equal to 2d _s in our model, i.e., 20% channel width or 0.4b.

Integrating the above algorithm, the model can detect incipient cutoffs. As described in earlier sections, the present model uses an irregular mesh based on triangular elements for the main channel and the floodplain; thus, it is crucial to keep track of the bank lines during meander evolution. Therefore, we define a binary variable to track the channel–floodplain regions. For the initial mesh, all channel nodes are given the value of 0, and all floodplain nodes are given the value of 1. The boundaries of “0 region” and “1 region” are two polylines representing the left and right bank lines (see Figure 3A). After several time steps, especially at the moment of cutoff, the “0 region” will be enlarged and forms a new left bank line or a right bank line (see Figure 3B). Next, the model will continuously run using the new configurations. Flow velocities in the newly-formed oxbow lake will reduce when time advances; thus, the “1 region” will gradually shrink and eventually form a new channel (see the process in Figures 3C–F).

Three categories of physical processes are involved: hydraulics, sediment transport/bed evolution, and migration of the river banks. Each of them has different timescales and for that reason, different time steps can be used for the numerical solution of the flow field, the bed evolution, and the computation of the migration of banks. The decoupling of the hydraulics and bed evolution is handled internally by SISYPHE. The time step is specified in term of multiples of the time used for the hydraulics (one, two or three orders of magnitude are typically utilized). The time step for the migration of the banks is managed by MEANDRE; the time step for the migration of the banks is set as several times the time step used to evolve the river bed. In cases where local processes are acting (e.g., bank collapse) at different rates than reach-averaged rates, a coupling between the different modules has to be accounted for.

The validation of the methodology formulated here was performed first using the results of the experiment of Han and Endreny (2014a); they utilized a physical model to simulate the development of a cutoff. The experiment was set in a mobile bed river table (3.66 m × 1.83 m × 0.2 m) where the water surface and riverbed were periodically monitored before, during and after the cutoff development. The water surface and riverbed configuration were measured at a sub-millimeter scale at various meander evolution stages: M3 was the initial stage; M4 was the pre-cutoff stage; M5 was at cutoff; M6 was 5 h after cutoff; all the mentioned stages are depicted in Figure 4A. The experiment used an 80 mL/s constant rate input of water into the river table’s upper reservoir while draining the lower reservoir from a fixed elevation orifice weir crest. The river had an average channel width and depth of 80 mm and 25 mm, respectively, a valley slope of 1.8%, a sediment size D₅₀ of 0.2 mm with 10% talc powder to create cohesivity, a mean velocity of 0.04 m/s, a Reynolds number of 3200, a Froude number of 0.05, an average boundary shear stress of 0.18 Pa and the bed material had a critical shear stress of 0.15 Pa. The river water surface and riverbed digital elevation model (DEM) data during the experiment were obtained from post-processing stereo-pair images (Han and Endreny, 2014b).

The cutoff modeling approach described herein was utilized to simulate this experiment and the results of flow velocity, bed shear stress, water depth and bed topography are shown in Figure 4B for the experimental meander stages M4, M5 and M6, which corresponded with the model Cycle 0, 8 and 24. A cycle is a single timestep where the river banks are adapted by the processes of bank erosion and accretion, and the mesh is generated again if necessary for a better representation of the river banks position. The model also simulated beyond M6 (Cycle 89), corresponding to a stage when the cutoff had transitioned to its lowest sinuosity. The water depth and velocity vectors at M5 (also shown at Figure 4) capture the development of the cutoff breaching the channel outer banks, and the changes in water depth between the stages M5 and M7 capture the transition from oxbow lake to paleo-channel. During cutoff, shear stresses increased in the cutoff region; these larger shear stresses then spatially expanded and cut a deeper, straighter channel through the meander neck, which resulted in their decrease.

The model simulated the post-cutoff sequence through Cycle 89 to track activity in the upstream erosional and downstream depositional processes, extending beyond M7, after the laboratory experiment terminated. This additional tracking allowed us to outline a conceptual model of the development of the longitudinal bed profile during and after the cutoff (Figure 5), which is based on the river geometry and other bend wavelengths and amplitudes found during cutoff evolution (Gutierrez and Abad, 2014; Gutierrez et al., 2014). The first stage of the conceptual model involves the conditions just before the cutoff, with the channel planform and longitudinal bed profile described in Figure 5A, where S _e denotes the equilibrium streamwise slope. The next stage is the development of the cutoff (Figure 5B), where the longitudinal bed profile is modified with a reach of slope, S _AB , where the shorter length of the new channel compared to the pre-cutoff conditions produced a larger slope. The next stage is indicated in Figures 5C, D; the channel tends to return to an equilibrium bed profile by triggering an erosional wave traveling upstream and a depositional wave traveling downstream, as indicated in Figure 5C, and depicted in the bed longitudinal profile; this process also promotes planform migration (Abad et al., 2012). A fundamental part of the conceptual model is the erosional and depositional waves (Figure 5C); the erosional wave moves upstream, incising the bed, and the depositional wave moves downstream, depositing the sediment eroded from the upstream area. This process helps the development of the oxbow lake with the definition of the new left bank line, as identified in Figure 5D. Moreover, the deposition of sediments downstream contributes to plugging the outlet of the oxbow lake.

The evolution of bed shear stress and bedload for the numerical modeling of Han and Endreny (2014a)’s cutoff experiment is shown in Figure 6. Bed shear stress increases after the cutoff, first locally near the neck, and upstream and downstream of it (see M5 and M6 in Figure 6A). After the meander loop is abandoned, the shear stress in the new shorter reach is larger than before the cutoff (compare M4 to M7 in Figure 6A). Note that after the cutoff, the downstream reach widens, because of the increased shear stress. Bedload upstream and downstream the cutoff increases to a peak (by cycle 30), and then gradually decreases (Figure 6B). The increment is larger downstream the cutoff. After the meander loop has been abandoned, and the reach of the cutoff has been straightened, bedload is larger than before the cutoff.

Following the model application in the laboratory scale, we conducted the simulation on a natural cutoff in the Ucayali River. In the satellite imagery provided by the Landsat Missions, a massive number of cutoffs, oxbow lakes, scroll bars, and paleo-channels can be observed along the Ucayali River (see the example reach illustrated in Figure 1). A cutoff event that occurred between 2014 and 2015 was considered for modeling with the approach presented in this study. The cutoff was located approximately 200 km downstream of Pucallpa City, eastern Peru. From the satellite imagery of 2006, 2010, 2013, 2014, 2015, 2016, and 2017, it was determined that the channel banks migrated at approximately 50–200 m/yr (see Figure 1). To set up the cutoff model based on the 2013 river geometry, several assumptions were made: 1) Since the goal was to reproduce the main physical processes of the cutoff (bank migration, riverbed erosion and sedimentation), the initial river geometry was considered of constant width of 800 m instead of the spatially varying river widths. The width was eventually self-adjusted during the simulation; 2) Due to the limitation of available data, the initial river bathymetry was obtained from a synthetic riverbed topography model proposed by Li and García (2021), which represents the morphodynamic bed equilibrium configuration around river bends. The computed bed topography was validated against the measured cross-sectional profiles. Similar to the channel width, the bed evolution will self-adjust the bathymetry during modeling; 3) The topography on the floodplain was not considered in this study. Hence, instead of applying local DEM data as the initial floodplain topography, the initial floodplain topography was modeled as an inclined plane with a constant valley slope; 4) A constant flow discharge of 12,000 m³/s (annual mean) was imposed as boundary condition at the inlet; 5) For the sediment boundary condition, a model of recirculating boundary condition was presented by Mendoza et al. (2017), where the volume of sediment measured at the outlet is fed back at the inlet. In such conditions, the sediment transport is at the dynamic equilibrium state, with oscillations caused by the development of bedforms. However, for the case of the Ucayali River, equilibrium conditions were set at the inlet sediment boundary condition, since after the neck cutoff, the sediment would have been fed at the inlet, something not necessarily true when only one cutoff is simulated.

Figure 7 shows the modeled migration of the banks between 2013 and 2017 in the analyzed reach; at Cycle 35, the cutoff was initiated. An abrupt drop in flow velocity is observed at the bend, that is, being abandoned. In Cycle 60 of the simulation, which corresponds chronologically to the satellite imagery of 2015, an oxbow lake was fully formed since the velocity is zero although the oxbow is still not infilled. There is a coincidence in the development of bars when comparing the modeling results; the bars marked with a to c indicated in the results of cycle 60 illustrated in Figure 7, are similar to the ones observed from the satellite imagery of 06/2015. During Cycle 100 the simulation developed a large bar just downstream of the cutoff point, which reconfigures the flow distribution and promotes bank erosion in front of such bar. This is also observed in the satellite imagery (image of 2017 in Figure 7). After the cutoff developed, the new shortened channel will progressively self-adjust to a new meander bend, whose channel sinuosity will gradually increase due to channel migration, same as the rest of the upstream and downstream bends (Figure 8). The simulation took 120 h for the Ucayali River in a high-performance cluster using 32 cores of computation. The channel (floodplain) was discretized using triangular elements with maximum size of 50 m (250 m) edge length.

The process of bed evolution also determines how the cutoff bend interacts temporarily with the flow before being completely abandoned. Such interaction is observed in the modeled flow velocity and bed elevation (Figure 7), where the velocity after the cutoff gradually diminished while the original bed was aggrading. The abandonment of the original bend occurs progressively. For instance, Turnipseed et al. (2021) showed the case of a neck cutoff in the White River, Arkansas, where during a bank full event the original meander still captured 51% of the flow 1 year after the cutoff.

The morphodynamic processes prior, during and after cutoff processes involve several complex patterns which can be addressed by the present modeling framework.

Paleochannels and river migration: As observed in Figure 1, during the lifetime of a river, both recent geomorphic features (e.g., scroll bars and oxbow lakes) and ancient geomorphic features (e.g., paleochannels and terraces) are found. The floodplain is usually characterized of fine sediment deposits and vertical layers of heterogenous soils (Abad et al., 2022), vegetation responding to river dynamics, differentiated elevations throughout the floodplain, and complex connectivity due to the existence of dry and wet paleochannels (Rojas et al., 2023), a required ingredient for preservation of the stratigraphic record in fluvial meander-belt deposits (Durkin et al., 2017).

Forced and migrating free bars: the coexistence of forced and free bars in meandering channels (Tubino and Seminara, 1990) is typical of meandering rivers. After cutoff, the channel develops a low sinuosity (less curvature forcing), and free migrating bars can interact with bars forced by channel enlargements (Luchi et al., 2010) as those produced by a neck cutoff (Figure 7). In addition the migration of both large scale (bars) and small scale (dunes) bedforms produces a modulation of the shear stresses exerted on the outer bank and, hence, on bank erosion (Abad et al., 2013).

Recently, detailed field measurements in cutoffs were performed along different rivers; White River, Arkansas (United States) (Richards et al., 2018; Richards and Konsoer, 2020; Turnipseed et al., 2021), Wabash River, Indiana and Illinois (United States) (Zinger et al., 2011), Huallaga River, Loreto (PERU) (Valverde et al., 2023). However, still a need for more detailed high-resolution hydrogeomorphic characterization of the different phases during pre-, cutoff, and post-cutoff conditions, including those of traveling sediment waves. In Amazonian rivers, cutoffs occur in a timescale dependent of the river characteristics, but the Ucayali River is one of most dynamic rivers where little hydrogeomorphic data is available (Abad et al., 2022). In order to improve river modeling at these large and very active systems, more detailed river data is needed, not just to understand the physical processes, but also to calibrate sediment transport equations, since, along these rivers, sediment fluxes are dependent on planform conditions (including influence of tributaries), flow variability and bedforms activity (Abad et al., 2022; Abad et al., 2023; Valverde et al., 2023).