Long-term hydro-sedimentary dynamics of the Ucayali River (Amazon Basin) revealed through combined observations, remote sensing, and SWAT-Amazon modelling

The Amazon basin is undergoing increasing environmental changes, potentially approaching a climatic tipping point in the coming decades. Understanding how these changes affect water and sediment fluxes is key for constraining large-scale biogeochemical cycles, yet conventional hydrological networks lack the spatial and temporal resolution required to accurately quantify hydro-sedimentary budgets.

To address this limitation, we develop an integrated, physically constrained framework combining long-term observations, remote sensing, and hydrological–hydraulic modelling (SWAT-Amazon) to quantify multi-decadal hydro-sedimentary budgets and investigate how floodplain inundation controls sediment dynamics in large Amazonian rivers. Focusing on the Ucayali River, a major foreland tributary of the Amazon, this study provides the first detailed, long-term hydro-sedimentary budgets for the Upper Amazon, distinguishing fine sediment fluxes from sand loads.

Results reveal a previously undocumented floodplain-controlled sand sedimentation process: during high waters, large floodplain water storage (up to 19.1 [15.3, 22.9] km³, ∼ 38 % of discharge) reduces main-channel transport capacity, capturing up to 14 % [10 %, 20 %] of the sand flux at peak discharge, while recycling during recession contributes 22 % of the total suspended load at the basin outlet. This dual control partially decouples sediment transport from water discharge. The Andean Ucayali exports 455 [410, 500] × 10⁶ t yr⁻¹ of suspended sediment (40 % sand), of which 36 % is trapped within the floodplain, predominantly as sand (65 % of total deposition). The river delivers 290 [235, 345] × 10⁶ t yr⁻¹ to the Amazon River (26 % sand), making it the dominant sediment source among the Andean foreland tributaries. Uncertainty analysis combining Sobol indices and GLUE simulations shows that, despite substantial equifinality among secondary floodplain parameters, sediment fluxes and associated trapping and recycling fractions remain stable across all behavioural simulations. Budget accuracy is therefore controlled by long-term, multi-variable, multi-source observations rather than by parameter calibration or model structure alone.

These findings demonstrate that floodplains control hydro-sedimentary fluxes in large river systems and act as dynamic regulators of sediment transport, storage, and recycling, with major implications for biogeochemical cycles.

1 Introduction

1.1 Global contribution of the Amazon Basin

The Amazon basin is a massive hotspot for water and matter inputs to the Ocean (Syvitski et al., 2005; Martinez et al., 2009; Moquet et al., 2016; Jouanno et al., 2021; Louchard et al., 2021, 2023) and plays a key role in global hydro-biogeochemical cycles (Gaillardet et al., 1999; Bouchez et al., 2012), capable to significantly impact oceanic biogeochemistry (Jouanno et al., 2021; Louchard et al., 2021). Long-term monitoring by the CZO (Critical Zone Observatory) HyBAm (Hydrology of the Amazon Basin) shows that the Amazon River annually discharges 6500 km³ of freshwater (∼ 20 %–25 % of the global total) (Callède et al., 2010), 1100×10⁶ t of suspended sediments (∼ 8 % of global riverine outputs) (Santini, 2020) and 272×10⁶ t of dissolved matter (Moquet et al., 2016) (∼ 7 % of the global flux). It also influences atmospheric circulation, contributing up to 15 % of global continental evapotranspiration (Salati, 1979; Soares-Filho et al., 2010; Satyamurty et al., 2013) and acts as both a carbon sink and a greenhouse gas source, contributing substantially to global cycles (Richey et al., 2002; Melack et al., 2004; Subramaniam et al., 2008; Ward et al., 2016; Pangala et al., 2017; Louchard et al., 2021). Nutrient-rich from the Andean Cordillera, the Amazon hosts 25 % of terrestrial species and the Earth's largest rainforest (e.g. Lesack, 1993; Malhi et al., 2008; Fan and Miguez-Macho, 2010).

1.2 Role of the floodplain dynamics

Amazonian floodplains act as dynamic reactors, playing a key role in global water and sediment fluxes. Lateral exchanges between the main channel and alluvial plains are of the same order of magnitude as the fluxes reaching the ocean (Meade et al., 1985; Mertes et al., 1996; Dunne et al., 1998) and dominate the annual floodplain water balance (Rudorff et al., 2014a, b). Around 30 % of peak discharge transits through the floodplain (Richey et al., 1989; Lininger and Latrubesse, 2016), with highly variable pathways and residence times (from seconds to months). The flooded area covers 8 %–10 % of the basin (5–6 ×10⁵ km²) (Fleischmann et al., 2022). Thus, flood dynamics regulate the storage and exchange of sediments, nutrients, organic matter, pollutants, and living organisms (Aufdenkampe et al., 2011; Lewin et al., 2017). Sediment residence time varies from brief periods to tens of thousands of years (Mertes et al., 1996; Allen, 2008), depending on the floodplain's geomorphology, influencing chemical maturation processes essential to global biogeochemical cycles, including CO₂ consumption by silicate weathering (Guyot et al., 2007; Bouchez et al., 2012). Flexural basins adjacent to the Eastern Cordillera trap 40 %–50 % of Andean sediment exports (Guyot, 1993; Baby and Guyot, 2009; Armijos et al., 2013; Santini et al., 2014; Vauchel et al., 2017; Santini, 2020). Further downstream, sediment balances tend toward equilibrium between deposition and resuspension during flood recession (Santini et al., 2014; Espinoza-Villar et al., 2017), though influenced by the basin's structural Arches. Downstream of the Amazon-Madeira confluence, floodplains have remained only partially filled since the last glacio-eustatic lowstand (∼ 125 m below present sea level) and act as fine-sediment sinks, where large, shallow lakes retain overbank floodwaters (Tricart, 1977; Fleming et al., 1998; Park and Latrubesse, 2017).

1.3 The impacts of global and local changes

The Amazon Basin is undergoing a dramatic transition (Walling, 2006; Malhi et al., 2008; Davidson et al., 2012), facing pressures from deforestation for agriculture and pasture, resources extraction, and construction of hydroelectric projects (Finer and Jenkins, 2012; Latrubesse et al., 2017; Timpe and Kaplan, 2017; Chaudhari and Pokhrel, 2022). In recent decades, the Amazon Basin has also been affected by global climate changes, experiencing more frequent extreme floods (e.g., in 2009, 2012, 2014, 2015) and severe droughts (e.g. in 2005, 2010, 2023, 2024), with an increase in the amplitude of the annual flood wave (Davidson et al., 2012; Espinoza et al., 2012, 2013; Marengo and Espinoza, 2016; Nobre et al., 2016; Towner et al., 2020). Maximum flooding extent along the central Amazon has expanded by 26 % (Fleischmann et al., 2023), mechanically impacting key processes such as CO₂ and CH₄ outgassing. This warming-induced hydrological cycle strengthening is projected to continue in the coming decades (Hirabayashi et al., 2013; Langerwisch et al., 2013; Sorribas et al., 2016; Alfieri et al., 2017) and the rainforest could reach a tipping point by the second half of the century, potentially converting to savanna, particularly in the eastern and southern regions, or persisting in a degraded state (Lovejoy and Nobre, 2019; McKay et al., 2022; Flores et al., 2024).

1.4 Monitoring challenges and integrated approach

The cascading effects of the ongoing transition in the Amazon remain uncertain, as monitoring is limited, particularly regarding sediment fluxes. Sparse measurements, due to high costs and logistical challenges, hinder the establishment of dense, long-term monitoring networks. Such networks are essential to constrain consistent upstream–downstream mass balances and to spatialize them at relevant scales. This is particularly important for identifying key processes, especially those linked to lateral exchange with the floodplain, which are still only roughly estimated. For instance, the CZO HyBAm covers just one gauging station per 160 000 km² on average. Additionally, sediment budgets in lowland sub-basins are challenging to estimate accurately, as their order of magnitude is comparable to the uncertainty in sediment load measurements (e.g. Xiaoqing, 2003; Horowitz et al., 2015; Vauchel et al., 2017; Gitto et al., 2017; Santini et al., 2019; Santini, 2020; Dramais, 2020). A significant portion of the suspended load (up to 70 %) consists of very fine sands (Santini et al., 2019; Martinelli, 2022), which are difficult to measure due to their sensitivity to hydrodynamic fluctuations and heterogeneous distribution within the cross section. Furthermore, the buffering effects of such a large basin (Walling, 2006) can mask the impacts on material transfer to the oceans, requiring long-term monitoring.

The scarcity and heterogeneity of observed data directly reduce the robustness and accuracy of hydrological and sediment transport models, limiting their ability to capture key processes and to reliably forecast responses to environmental changes. In response, spatial hydrology has increasingly complemented in situ observations in the Amazon, with satellite data, particularly space altimetry and water color imaging, playing a key role in monitoring water levels and sediment concentrations (Calmant et al., 2009; Martinez et al., 2009, 2015; Espinoza-Villar et al., 2013, 2017; Park and Latrubesse, 2014). On the other hand, hydrological models have also addressed flooding and backwater effects (Yamazaki et al., 2011; Paiva et al., 2013; Pontes et al., 2017; Siqueira et al., 2018; Santini, 2020; Guilhen et al., 2022). Recently, the question of the sediment routing into semi-distributed models has been explored (Fagundes et al., 2021, 2023; Santini, 2020).

However, to date, no study has yet combined remote sensing with modelling to investigate sediment dynamics in detail in the Amazon. Building on long-term CZO HyBAm conventional observations, this study introduces an integrated framework coupling field calibration campaigns, satellite remote sensing, and hydraulic–hydrological modelling to derive process-based hydro-sedimentary budgets for a major foreland tributary of the Upper Amazon: the Ucayali River. Satellite altimetry and satellite-derived fine sediment estimates are used to constrain a hydrological–hydraulic modelling scheme (SWAT-Amazon) simulating water and sand fluxes, with floodplains represented using a simplified reservoir approach.

The objectives are to (i) quantify long-term hydro-sedimentary budgets at sub-basin scale, (ii) investigate the impact of floodplain inundation on water and sediment fluxes, thereby identifying key processes controlling sediment transport, storage, and recycling that remain poorly constrained by conventional river monitoring approaches, and (iii) distinguish the respective contributions of fine sediment fluxes, associated with organic matter and pollutant transfer, and sand loads related to river dynamics.

Unlike previous basin-scale sediment studies in the Amazon, which largely relied on sparse gauging networks or large-scale modelling approaches, this framework provides process-based hydro-sedimentary budgets through multi-source data integration. Given the key role of sediment dynamics in biogeochemical cycles, it also contributes to improving the understanding of the Amazon's role in global material fluxes and of the potential impacts of environmental changes on its hydrology and sediment transport.

1.5 Case study: the Ucayali Basin

Given the continental scale of the Amazon Basin, this study focuses on the Ucayali River, a major foreland tributary draining 350 000 km² (49 % Andes, 51 % plains). Only two HyBAm gauging stations monitor hydro-sedimentary fluxes in this basin, with incomplete records for the study period (1983–2019): one upstream of the lowlands, the other at the basin outlet. Limited water levels, with some unreliable records, along with a few discharge measurements, are also available from other lowland stations. These constraints make the Ucayali a relevant test site for building the proposed integrated approach, before extending it to other Amazonian sub-basins.

2 Integrative strategy

The tailored integrative strategy (Fig. 1) for improving water and sediment balances relies on a combination of three station types in the plain: (i) “low-data” (or poorly monitored) conventional stations, characterized by incomplete and/or inconsistent datasets; (ii) “virtual” stations established at locations where satellite altimetry ground tracks intersect the river mainstem, in order to enhance the spatial density of the monitoring network through the integration of remote sensing and modelling; and (iii) “super” stations with long-term, high-quality datasets, which serve as benchmarks for calibrating and validating the integrated approach. This strategy is further supported by dedicated calibration field campaigns. At all stations, water and sediment fluxes are estimated by integrating remote sensing products and hydrological modelling outputs. Water discharges are simulated using a modified version of the Soil and Water Assessment Tool (SWAT) model (Arnold et al., 1998), to account for Amazon flood wave dynamics and attenuation during flooding. It is assumed that the river transports two particle groups (Santini et al., 2019): fine sediments and sands (Fig. 1b). Fine sediments (mean diameter d_f≅ 10–20 µm), primarily silts with small clay aggregates, behave similarly to passive scalars and their fluxes are not modelled with transport capacity equations. Instead, fine sediment concentrations at the water surface are derived from satellite images (Espinoza-Villar et al., 2012, 2013, 2017; Martinez et al., 2009, 2015), using an inversion model calibrated with in situ data. Whereas suspended sands (mean diameter d_s≅ 80–120 µm) transported in graded suspension, are invisible to satellite spectroradiometers due to Mie scattering (Pinet, 2017). These ranges of mean diameter should be understood as representative values of the dominant PSD modes, rather than strict grain size distribution bounds. Intermediate particle sizes (e.g. coarse silts in the 20–63 µm range) are included within the fine sediment fraction defined by the 63 µm threshold used in the HyBAm monitoring protocol (see Sect. 3.1). Depending on their physical properties, coarse silts may either behave as wash load (e.g. aggregates) or exhibit transport dynamics closer to those of very fine sands when non-cohesive. This supports a first-order bimodal representation of suspended sediment transport at the basin scale, consistent with PSD measurements and vertical grain-size profiles presented in the Supplement (Fig. S1).

Figure 1 General schematic overview of the proposed methodology. Panels (a–c) illustrate the integrative approach. (a) Types of stations. (b) Typical bimodal particle size distribution (PSD) in the large Amazonian rivers, identifying two main size groups: 1 – fine sediments that can be monitored by satellite but not modelled; 2 – fine sands in graded suspension, invisible to satellites but whose transport capacity can be modelled. (c) Integrated approach combining remote sensing, modelling, and calibration campaigns. (d) SWAT-Amazon, a tailored version of the SWAT model for simulating water and sand fluxes. This modelling framework consists of a Fortran-based executable (SWAT-Amazon.exe) and an R notebook (Run-SWAT-Amazon.Rmd) used for model runs, simulation analysis, interactive visualization, sensitivity analysis and calibration with the SWATrunR package (Schürz, 2019).

Moreover, according to Santini et al. (2019) and Martinelli (2022), observed Rouse numbers (Rouse, 1937) are between 0.2 and 0.8 for this sand fraction, inducing concentrations near the surface. Therefore, sand loads are modelled using sediment transport equations in a new routing module developed in the SWAT model, referred to as SWAT-Amazon (Fig. 1d).

3 Dataset

3.1 Conventional data

This study relies on long-term hydro-sedimentary flux data from the CZO HyBAm (Guyot et al., 2007; Santini, 2020). In the Ucayali basin, IRD (Institut de Recherche pour le Développement) and SENAMHI (Servicio Nacional de Meteorología e Hidrología) have been operating two HyBAm gauging stations, Lagarto (HyBAm code: 10073500) and Requena (10074800) since 2001 (Fig. 2), carrying out 82 field campaigns to establish rating curves. Additional sediment monitoring was carried out at Puerto Inca between 2012 and 2016, at a conventional SENAMHI station. Water level–discharge relationships were also established at Puerto Inca (10073750) and Pucallpa (10074000), where water levels are monitored by port authorities. Field measurement protocols are described in Sect. S2 in the Supplement. Observed sand fluxes, empirically derived from gauging and surface concentration (S3 in the Supplement), carry ±30 % uncertainties, affecting simulations statistics.

Figure 2 Ucayali's stack (a) Ucayali Basin sub-basins and station locations. Super stations (CZO HyBAm) (blue-filled circles, white outline), virtual stations (green-filled circles, blue outline), and low-data Stations with sparse/inconsistent data (blue-filled circles, green outline). The text box details observation periods for water level (h) discharge (Q), and suspended sediment concentration (C). Main Andean tributaries on the left bank: Cushabatay (3); Pisqui (9); Aguaytia (8); Pachitea (10, 20, 6). (b) Geomorphological domains and outcrop distribution. Box 1: Narrow sedimentary basin controlled by Shira Mountains and Fitzcarrald Arch uplift; Box 2: Piggyback Basin backward the Moa Divisor Thrust Fault; Box 3: Marañon Foredeep. Wetlands extent based on Hess et al. (2015). (c) Mean annual precipitation map (PISCO dataset, 1983–2019).

3.2 Altimetric data and definition of virtual stations

Two virtual stations, JA204-S3B310-R22 (reach 22 on Fig. 2) and S3B331-R5 (reach 5), were defined based on the intersection of satellite altimetry ground tracks (Jason, Envisat, Sentinel) with the mainstem of the Ucayali River. These stations provided satellite-derived water level time series used to calibrate the hydrological model. In addition, satellite altimetry was employed to correct water level records at Requena, Pucallpa and Contamana (10074500), the latter being a rarely visited SENAMHI station without any flow measurements. Altimetry data processing was carried out using the open-access VALS (Virtual ALtimetric Stations) software. By summing the virtual stations, the low-data stations (Contamana and Pucallpa), and the long-term CZO HyBAm stations (Lagarto and Requena), the Ucayali sedimentary basin was subdivided into five distinct compartments to establish hydro-sedimentary budgets with the integrated approach.

3.3 Fine sediments monitoring with remote sensing data

3.3.1 Retrieving time series of remote-sensed reflectance data

Given the required revisit frequency and the Ucayali River's width in the plains (500–1000 m), moderate-resolution satellite imagery from MODIS (MODerate Resolution Imaging Spectroradiometer, 250 × 250 m, 1999–present, 1–2 d) and VIIRS (Visible Infrared Imaging Radiometer Suite, 375 × 375 m, 2012–present, 0.5 d) was used to generate time series of surface water reflectance (Fig. 1c). Reflectance values in the red and near-infrared (NIR) bands were extracted pixel by pixel from satellite images using the free software GetMODIS and MOD3R, developed by the CZO HyBAm, that was tested and validated in various previous studies (e.g. Espinoza-Villar et al., 2017; Vauchel et al., 2017). Water masks were applied to the Ucayali River's main course near virtual stations. To ensure 50–100 pixels per mask, large river stems were covered, with masks redrawn every 2–3 years due to river mobility. Collected scenes comprise images spanning 8 d periods, selecting pixels with the lowest cloud cover and smallest satellite-viewing nadir angle.

3.3.2 Conversion of remotely sensed reflectance to fine sediment concentration

Two radiometric campaigns were conducted in the Ucayali Basin: the first in November 2011 at Requena (Espinoza-Villar et al., 2012) and the second in February 2017 at Lagarto, Puerto Inca, and Pucallpa, spanning three weeks (Santini, 2020). A total of 42 surface water samples were collected to determine total, fine, and sand concentrations. Simultaneously, hyperspectral field radiometers (TriOS) were deployed following the experimental setup of Mobley (1999), as adapted by Martinez et al. (2015) for the Amazon Basin. High-frequency (1 Hz) hyperspectral measurements of surface water reflectance were obtained at sampling locations. Relying on this dataset, a unique model for all the Ucayali Basin was fitted between fine sediment concentration at the water surface and the ratio of radiometer reflectance in the NIR (841–876 nm, according to the satellite sensor bands) and red bands (620–670 nm) (see Sect. 5.4). This single relationship is considered applicable along the mainstem of the Ucayali River in the lowland plain because surface reflectance is largely controlled by fine silts of Andean origin, whose grain-size characteristics and optical properties are relatively homogeneous along the river continuum (Martinez et al., 2015; Santini, 2020). In addition, the use of spectral band ratios (NIR $/$ RED), commonly applied to reduce the influence of these potential variations in optical conditions (Doxaran et al., 2002; Martinez et al., 2015; Pinet, 2017), supports the applicability of this relationship at the basin scale considered in this study.

3.3.3 From surface to mean concentration of fine sediments

Due to the considerable depth of Amazonian rivers and the vertical sediment concentration gradient near the surface, the ratio α_f, relating the channel mean concentration to the surface index concentration retrieved by satellite, ranges from 1 to 1.8 according to the CZO HyBAm database (1.1 to 1.2 in the Ucayali). These values are derived from field measurements of vertical suspended sediment concentration profiles collected at several stations of the HyBAm observatory across the Amazon Basin, including the Ucayali River (e.g. Santini et al., 2019; Santini, 2020). To estimate α_f along the river network, the models proposed by Santini et al. (2019) were applied. These models use Rouse-type formulations constrained by observed concentration profiles to describe the vertical distribution of suspended sediments as a function of hydraulic conditions. They were parameterized using hydraulic variables simulated by SWAT-Amazon.

3.4 Input data for modelling

The study utilizes the Peruvian Interpolated Data of SENAMHI's Climatological Observations (PISCO) (Aybar et al., 2020; Llauca et al., 2021) to support the development of an operational model in collaboration with the SENAMHI. Potential Evapotranspiration (PET) was estimated using the Hargreaves and Samani (1985) method, to take advantage of PISCO's temperature data. Land use data was obtained from the Peruvian Ministry of Environment (https://geoservidor.minam.gob.pe, last access: 6 June 2025), while soil information was sourced from the Harmonized World Soil Database (https://fao.org/soils-portal, last access: 6 June 2025). The topography layer was derived from the Multi-Error-Removed Improved-Terrain Digital Elevation Model (MERIT DEM) (Yamazaki et al., 2017), resampled from 90 to 300 m for computational efficiency.

In SWAT-Amazon, water and sand fluxes can be forced at any sub-basin via input files. However, no external forcing was applied in this study. For sand fluxes, we assumed that Andean inflows (basins 3, 8, 9 in Fig. 2) were governed solely by transport capacity within the sand routing module. This assumption is supported by the observed relationship between sand flux and water discharge at Lagarto and Puerto Inca, which indicates, to a first approximation, sediment availability throughout the hydrological cycle in the Andean sub-basins. Lateral contributions from plain tributaries (basins 1, 2, 4, 7) were considered negligible and were likewise represented as transport-capacity limited in the simulations.

4 Tailoring SWAT for water and sediment flux simulation

SWAT, a semi-distributed model with physical and conceptual equations, was chosen for its proven robustness in simulating hydrological processes in large basins at daily time step. Its open-source Fortran code and extensive user community provide numerous complementary modules and tools. However, SWAT has limitations in modelling water and sediment routing in large rivers with diffusive flood waves and extensive floodplains. It lacks realistic hydraulic connectivity between floodplains and the main channel, preventing accurate simulation of the relationships between water levels, velocities, and discharges, which are keys for sediment transport. To address this, a major code modification is introduced below.

4.1 New water routing modules

The main channel's trapezoidal cross-section in SWAT was replaced with a rectangular one for consistency with the hydraulic equations used. Floodplains were represented using simplified rectangular or triangular cross-sections. Alternative geometries were tested during model development (Santini, 2020), and the selected configuration was retained as the most robust and parsimonious for large low-slope rivers. In this framework, channel geometry parameters represent effective hydraulic conditions averaged at the reach scale rather than the instantaneous geometry of a single cross-section. Although large Amazonian rivers are morphodynamically active and may exhibit lateral migration rates of several tens of meters per year, long-term observations from the HyBAm network indicate that stage-discharge relationships remain remarkably stable over time at the mainstem stations. This suggests that flow conditions are primarily controlled by reach-scale channel geometry over the decadal time scales considered here.

4.1.1 Water level calculation and state variables

The mean water level h (m) is derived from the water volume V (m³) stored in the reach i at time step t (beginning of the simulation day). As long as h≤h_f, where h_f (m) is the floodplain activation threshold, h is computed as:

$$\begin{array}{}\text{(1)}& h_i^t={\displaystyle \frac{V_i^t}{\mathrm{\Delta }{x}_i{B}_i}},\end{array}$$

where B (m) is the main channel width and Δx (m) is the reach length. When h≥h_f, the floodplain is activated, distributing V between the main channel and floodplain. If the floodplain cross-section is rectangular:

$$\begin{array}{}\text{(2)}& h_i^t={\displaystyle \frac{V_i^t-\left(\mathrm{\Delta }{x}_i{B}_i{h}_{\mathrm{f},i}\right)}{\mathrm{\Delta }{x}_i\left(B_i+W_{\mathrm{fp},i}\right)}}+h_{\mathrm{f},i},\end{array}$$

where W_fp=k_fpB (m) is the floodplain width and k_fp (–) a coefficient to be calibrated (the SWAT-Amazon parameters are given in Table 1). For a triangular floodplain cross-section, h depends on θ_fp (rad), the riverward slope angle (see S4 in the Supplement). Finally, state variables such as wetted area A_h (m²), wetted perimeter P_h (m), and hydraulic radius R_h (m) are all derived from h.

Table 1 Main parameters in the new routing modules of SWAT-Amazon.

4.1.2 Dynamic and process equations

In large Amazonian rivers, flow variations over time and space are minimal, leading to a subcritical hydraulic regime, that can be modelled using the 1D Barré de Saint-Venant equations (Moussa and Bocquillon, 2009). Given the (very) gradual flow variations (Trigg et al., 2009), the convective and local acceleration terms are negligible, making the diffusive flood wave approximation suitable. When water-surface slope effects are also negligible, the pressure gradient is eliminated, allowing the use of the kinematic wave equation. SWAT-Amazon enables reach-specific selection between kinematic wave (S_f=S_b), suitable for steep Andean reaches, and diffusive wave approximation ($S_{\mathrm{f}}=S_{\mathrm{b}}+S_{\mathrm{w}}$), preferred for low-slope floodplain reaches where backwater effects may occur (e.g. Yamazaki et al., 2011), with S_f (m m⁻¹) the energy gradient (or friction slope), S_b (m m⁻¹) the bed slope and S_w (m m⁻¹) the water-surface slope. When the diffusive wave model is used, S_w is assessed as follows:

$$\begin{array}{}\text{(3)}& S_{\mathrm{w},i}^t={\displaystyle \frac{\left(h_{i+\mathrm{1}}^t-h_i^t\right)}{\frac{\mathrm{1}}{\mathrm{2}}\left(\mathrm{\Delta }{x}_i+\mathrm{\Delta }{x}_{i+\mathrm{1}}\right)}}.\end{array}$$

The reach-averaged velocity u (m s⁻¹) and discharge Q (m³ s⁻¹) are then calculated using the Gauckler-Manning-Strickler (GMS) friction equation (Hager, 2005).

$$\begin{array}{}\text{(4)}& Q_i^t=A_{\mathrm{h},i}^t{u}_i^t=A_{\mathrm{h},i}^t{\displaystyle \frac{\sqrt{S_{\mathrm{f},i}^t}}{n_{\mathrm{c},i}^t}}{\left(R_{\mathrm{h},i}^t\right)}^{\mathrm{2}/\mathrm{3}}.\end{array}$$

In the present study, the diffusive wave formulation was applied throughout the river network in order to account for backwater effects in the low-slope reaches of the Ucayali River.

4.1.3 Continuity equation and water storage in the reach

At the end of the calculation time step (t+Δt), V in reach i is updated as:

$$\begin{array}{}\text{(5)}& V_i^{t+\mathrm{\Delta }t}=V_i^t+\left(Q_{i-\mathrm{1}}^t-Q_i^t\right)\mathrm{\Delta }t+{V_{{\mathrm{R}}_i}}^t-{V_{{\mathrm{E}}_i}}^t-{V_{{\mathrm{T}}_i}}^t,\end{array}$$

where V_E (m³) is the volume lost to evaporation, V_T (m³) is the volume infiltrated into the unsaturated water table and V_R (m³) is the runoff (surface, subsurface, and baseflow) reaching the river. The computation of V_E, V_T and V_R follows the standard SWAT model. The updated volume $V_i^{t+\mathrm{\Delta }t}$ is then used to determine $h_i^{t+\mathrm{\Delta }t}$ for the next simulation step.

Although SWAT simulations are reported at a daily time step, river routing is internally solved at a sub-daily time step dynamically determined by the Courant–Friedrichs–Lewy (CFL) stability condition. The routing module automatically adjusts the internal time step to satisfy the CFL criterion and ensure numerical stability (Bates et al., 2010). Daily outputs therefore correspond to the model state at the end of each simulation day.

4.2 New module for sand sediment routing.

4.2.1 Sand load and concentration in the reach

At the beginning of the simulation day, the suspended sand concentration C_s (t m⁻³) in the reach i is:

$$\begin{array}{}\text{(6)}& C_{\mathrm{s},i}^t={\displaystyle \frac{V_{\mathrm{s},i}^t}{\frac{Q_i^t}{u_i^t}\mathrm{\Delta }{x}_i}},\end{array}$$

where V_s (t) is the sand volume stored in the reach, in the main channel only. The daily suspended sand load Q_s (t d⁻¹), taking Δt= 86 400 s, is then:

$$\begin{array}{}\text{(7)}& Q_{\mathrm{s},i}^t=Q_i^t{C}_{\mathrm{s},i}^t\mathrm{\Delta }t,\end{array}$$

4.2.2 Transport capacity evaluation

Selecting appropriate transport capacity equations for deep, low-gradient rivers is crucial, as most were derived from laboratory studies under opposite conditions (steep slopes, shallow water, uniform flow). The physically based Camenen and Larson (2005, 2008) models for non-cohesive sands were chosen for their calibration with extensive global datasets and proven applicability in large tropical rivers (Camenen et al., 2014). In these models, the transport capacity $Q_{\mathrm{s}}^{\ast }$ (t d⁻¹) for suspended sands is evaluated as a function of the Rouse number P_s (–), which defines the concentration profile exponential shape, and a near-bed reference concentration $C_{\mathrm{b}}^{\ast }$ (m³ m⁻³), which determines its magnitude:

$$\begin{array}{}\text{(8)}& {Q_{\mathrm{s}}^{\ast }}_i^t={C_{\mathrm{b}}^{\ast }}_i^t\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{6}{P}_{{\mathrm{s}}_i}^t}}\left(\mathrm{1}-\exp \left(-\mathrm{6}{P}_{{\mathrm{s}}_i}^t\right)\right)\right]Q_i^{t}s\mathrm{\Delta }t,\end{array}$$

where s is the relative sand density. The Rouse number summarizes the equilibrium between grain settling velocity w_s (m s⁻¹) and turbulence-induced lift, related to the shear velocity u_∗ (m s⁻¹), weighted by the sediment-to-eddy diffusivity ratio β_s (–):

$$\begin{array}{}\text{(9)}& {P_{{\mathrm{s}}_i}}^t={\displaystyle \frac{w_{\mathrm{s},i}}{{\mathit{\beta }}_{\mathrm{s},i}^t\mathit{\kappa }{u}_{{\ast }_i}^t}},\end{array}$$

where κ is the Von Kármán constant. The Soulsby (1997) law is used for estimating the sand grain settling velocity, involving the grain size d_s (m) of the suspended sands. The shear velocity is calculated using the depth-slope product:

$$\begin{array}{}\text{(10)}& {u_{{\ast }_i}}^t=\sqrt{g{h}_i^t\left({S_{{\mathrm{b}}_i}}^t+{S_{{\mathrm{w}}_i}}^t\right)},\end{array}$$

where g (m s⁻¹) is the gravitational acceleration. The diffusivity ratio is either assigned a fixed value for each reach or computed dynamically using the Santini et al. (2019) model:

$$\begin{array}{}\text{(11)}& {\mathit{\beta }}_{\mathrm{s},i}^t=\mathrm{3.1}\exp \left[-\mathrm{0.19}\times {\mathrm{10}}^{-\mathrm{3}}{\displaystyle \frac{u_{{\ast }_i}^t}{w_{\mathrm{s},i}}}{\left({\displaystyle \frac{h_i^t}{d_{\mathrm{s},i}}}\right)}^{\mathrm{0.6}}\right]+\mathrm{0.16},\end{array}$$

The bottom reference concentration $C_{\mathrm{b}}^{\ast }$ is given by Camenen and Larson (2005):

$$\begin{array}{}\text{(12)}& C_{{\mathrm{b}}_i}^{\ast t}={\displaystyle \frac{\mathrm{0.0015}{{\mathit{\theta }}^{\prime }}_i^t}{\exp \left(\mathrm{0.2}{d}_{\mathrm{b},i}+\mathrm{4.5}\frac{{\mathit{\theta }}_{\mathrm{cr},i}}{{{\mathit{\theta }}^{\prime }}_i^t}\right)}},\end{array}$$

where ${\mathit{\theta }}^{{}^{\prime }}$ is the dimensionless grain-related bed shear stress and θ_cr (–) the critical Shields parameter for the inception of transport (Camenen et al., 2014), which can be estimated from the Yalin-Shields curve as a function of the riverbed dimensionless mean diameter d_b∗:

$$\begin{array}{}\text{(13)}& {\mathit{\theta }}_{\mathrm{cr},i}={\displaystyle \frac{\mathrm{0.25}}{d_{\mathrm{b}\ast ,i}}}+\mathrm{0.055}\left(\mathrm{1}-\exp \left(-\mathrm{0.02}{d}_{\mathrm{b}\ast ,i}\right)\right).\end{array}$$

The parameter ${\mathit{\theta }}^{{}^{\prime }}$ is calculated as following:

$$\begin{array}{}\text{(14)}& {{\mathit{\theta }}^{\prime }}_i^t={\displaystyle \frac{\frac{\mathrm{1}}{\mathrm{2}}{f^{\prime }}_i^t{\left(u_i^t\right)}^{\mathrm{2}}}{g\left(s-\mathrm{1}\right)d_{\mathrm{b},i}}},\end{array}$$

where f^′ (–) is the Darcy-Weisbach skin roughness factor, derived from the logarithmic velocity law:

$$\begin{array}{}\text{(15)}& {f^{\prime }}_i^t=\mathrm{2}{\left({\displaystyle \frac{\mathit{\kappa }}{\ln \left(\frac{\mathrm{30}}{k_{\mathrm{s},i}^{\prime }}{h}_i^t\right)-\mathrm{1}}}\right)}^{\mathrm{2}},\end{array}$$

with κ is the Von Kármán constant and the height $k_{\mathrm{s}}^{\prime }$ (m) is the hydraulic skin roughness of Nikuradse, which can be expressed as a function of d_b and a coefficient σ=2.5 (Engelund and Hansen, 1967; Bartholdy et al., 2010):

$$\begin{array}{}\text{(16)}& k_{\mathrm{s},i}^{\prime }={\mathit{\sigma }}_i{d}_{\mathrm{b},i}.\end{array}$$

4.2.3 Sand load adjustment based on transport capacity

The difference $\mathrm{\Delta }{Q}_{\mathrm{s}}^{\ast }$ (t d⁻¹) between the transport capacity $Q_{\mathrm{s}}^{\ast }$ (t d⁻¹) and the sand load Q_s (t d⁻¹) is then calculated at time t:

$$\begin{array}{}\text{(17)}& \mathrm{\Delta }{Q}_{{\mathrm{s}}_i}^{\ast t}=Q_{{\mathrm{s}}_i}^{\ast t}-Q_{{\mathrm{s}}_i}^t.\end{array}$$

If $\mathrm{\Delta }{Q}_{\mathrm{s}}^{\ast }>\mathrm{0}$, there is excess transport capacity, allowing for riverbed erosion. The eroded mass E_bed (t d⁻¹) is defined as:

$$\begin{array}{}\text{(18)}& {E_{\mathrm{bed}i}}^t=K_{\mathrm{bed}i}\mathrm{\Delta }{Q}_{{\mathrm{s}}_i}^{\ast t},\end{array}$$

where K_bed (–) is a coefficient ($\mathrm{0}\le K_{\mathrm{bed}}\le \mathrm{1}$) representing the susceptibility of the channel bed to erosion when the simulated sand transport capacity exceeds the available sand load. This parameter governs the potential entrainment of bed material into suspension and thus the possible contribution of riverbed erosion to the simulated suspended sand flux.

The sand flux is then updated:

$$\begin{array}{}\text{(19)}& {Q_{{\mathrm{s}}_i}}^t={Q_{{\mathrm{s}}_i}}^t+{E_{{\mathrm{bed}}_i}}^t.\end{array}$$

Conversely, if $\mathrm{\Delta }{Q}_{\mathrm{s}}^{\ast }\le \mathrm{0}$, the sand load exceeds transport capacity, and the sand load is set to the transport capacity:

$$\begin{array}{}\text{(20)}& {Q_{{\mathrm{s}}_i}}^t=Q_{{\mathrm{s}}_i}^{\ast t}.\end{array}$$

4.2.4 Sand budget at reach scale

Drawing on the mass balance proposed by Dunne et al. (1998), an erosion term, E_bk (t d⁻¹), is introduced to account for both floodplain channel inputs and bank erosion. These processes primarily occur at point bars on the inner bends of meanders, where floodplain inflows, with lower sediment concentrations than the river's transport capacity, enhance erosion and resuspension. Riverbed erosion, E_bed (Eq. 18) and two deposition terms (D_ovbk, D_lat), all expressed in (t d⁻¹), are also considered. Riverbed erosion is therefore explicitly represented in the model formulation so that its potential contribution to the suspended sand flux can be evaluated during calibration. Sand deposition on bars in low-velocity zones of the main channel is already accounted for in the sand load adjustment (Eq. 20). Deposition in floodplain channels and levee depressions when active (i.e. when h>h_f) is neglected for sand particles, as the high flow resistance caused by vegetation in these areas is expected to result in complete sedimentation at their inlets. This process is therefore implicitly included in the overbank deposition term, D_ovbk.

The term E_bk is activated only when the daily water volume ΔV_fp (m³ d⁻¹) exchanged between the main channel and the floodplain is negative, meaning floodplain waters contribute to the main channel. Thus, E_bk was defined as function of ΔV_fp and C_bk (t m⁻³), the concentration of these banks and bars inputs, considered as a constant to be calibrated:

$$\begin{array}{}\text{(21)}& {E_{{\mathrm{bk}}_i}}^t=-\mathrm{\Delta }{V}_{\mathrm{fp},i}{C_{{\mathrm{bk}}_i}}^t.\end{array}$$

Thus, E_bk is neglected when ΔV_fp≥0 because the volume of water that could flow back from the floodplain during the rising stages is low compared to the water discharge in the main channel, contrary to the flood recession phase. In addition, the term E_bed can compensate for this if necessary, when the transport capacity is in excess.

The daily sand mass D_ovbk is defined as a function of ΔV_fp:

$$\begin{array}{}\text{(22)}& {D_{{\mathrm{ovbk}}_i}}^t=\mathrm{\Delta }{V}_{\mathrm{fp},i}{C}_{{\mathrm{s}}_i}^t\left({z_{{\mathrm{surf}}_i}}^t\right),\end{array}$$

where C_s(z_surf) (t m⁻³) is the sand concentration in the upper flow layer, estimated at $z_{\mathrm{surf}}\cong \left(h-h_{\mathrm{f}}\right)/\mathrm{2}$ (m). C_s(z_surf) is derived from the mean concentration C_s in the reach:

$$\begin{array}{}\text{(23)}& {C_{{\mathrm{s}}_i}}^t\left({z_{{\mathrm{surf}}_i}}^t\right)={\displaystyle \frac{{C_{{\mathrm{s}}_i}}^t}{{\mathit{\alpha }}_{\mathrm{s}}}}.\end{array}$$

The ratio ${\mathit{\alpha }}_{\mathrm{s}}=C_{\mathrm{s}}\left(z_{\mathrm{surf}}\right)/C_{\mathrm{s}}$ is estimated with the Santini et al. (2019) model:

$$\begin{array}{}\text{(24)}& {{\mathit{\alpha }}_{{\mathrm{s}}_i}}^t={\displaystyle \frac{\mathrm{1}}{\mathrm{6}{P}_{\mathrm{s},{\mathrm{cf}}_i}^t}}\exp \left(\mathrm{6}{P_{{\mathrm{s}}_i}}^t{\displaystyle \frac{{z_{{\mathrm{surf}}_i}}^t}{h_i^t}}\right)\left(\mathrm{1}-\exp \left(-\mathrm{6}{P_{{\mathrm{s}}_i}}^t\right)\right).\end{array}$$

4.2.5 Continuity equation

At the end of the calculation time step, the sand volume V_s stored in the main channel of the reach i is updated as:

$$\begin{array}{}\text{(25)}& V_{\mathrm{s},i}^{t+\mathrm{\Delta }t}=V_{\mathrm{s},i}^t+\left(Q_{\mathrm{s},i-\mathrm{1}}^t-Q_{\mathrm{s},i}^t\right)\mathrm{\Delta }t-D_{\mathrm{ovbk},i}^t+E_{\mathrm{bk},i}^t,\end{array}$$

As for water routing, sand routing is internally solved at a sub-daily time step to satisfy the CFL stability condition. The routed sediment inflow provided by the SWAT-Amazon reach network represents the sediment supply entering each reach and is incorporated into the reach-scale sand mass balance described in Eq. (25).

4.3 Additional flow resistances

4.3.1 Impact of floodplain activation on flow velocity and transport capacity

When the floodplain becomes active, differences in depth and roughness between the main channel and floodplain develop a shear interface between the two flow zones, associated with Kelvin-Helmholtz instabilities, transferring horizontal momentum from the main channel to the floodplain (e.g. Sellin, 1964; Nicollet and Uan, 1979; Ervine and Baird, 1982; Knight, 1989; Knight and Shiono, 1996; Smart, 1992; Loveless et al., 2000; Yen, 2002; Uijttewaal, 2014; Atabay and Knight, 2018; Proust and Nikora, 2020). Sediment-laden water flowing through floodplain channels (Lewin et al., 2017) also transfer large amounts of momentum to the plain and reduces the kinetic energy of the main flow, as does the attenuation of the water surface slope during flooding, which tends toward the valley slope. Moreover, the waters that travel for a short time through the floodplain before returning to the main channel also contribute to reduce the flow velocity. These combined effects significantly reduce flow velocity and, more drastically, transport capacity in the main channel. They change the spatial distribution of velocities and shear stress in the main channel cross-section, especially near the banks and bars, where sediment stocks can be available. To account for this, a flow resistance correction factor ζ_n was defined as:

$$\begin{array}{}\text{(26)}& {\mathit{\zeta }}_n={\displaystyle \frac{u_{\mathrm{cf}}}{u_{\mathrm{c}}}}={\displaystyle \frac{n_{\mathrm{c}}}{n_{\mathrm{cf}}}},\end{array}$$

where the subscripts “c” and “cf” denote in-bank flow configuration (without floodplain) and flow with an active floodplain, respectively, at the same water level h>h_f. To evaluate ζ_n, the Nicollet and Uan (1979) or Smart (1992) equations can be used. However, both formulations only consider the shear layer interface between the main channel and floodplain. Furthermore, the Smart equation is not suitable for large rivers, and the Nicollet et Uan equation requires an estimate of the floodplain's Manning coefficient. Although the latter was implemented in the new water routing module, a simpler approach was preferred. Therefore, a relationship between ζ_n and the relative height $Y=\left(h-h_{\mathrm{f}}\right)/h$ (–), from which the water exchanges between the main channel and the floodplain begin to affect the flow velocity, was defined:

$$\begin{array}{}\text{(27)}& {{\mathit{\zeta }}_{ni}}^t={\displaystyle \frac{\mathrm{1}}{\left(\mathrm{1}+Y_i^t{C}_{n\mathrm{fp},i}\right)}},\end{array}$$

with C_nfp (–) a coefficient to be calibrated, superior to zero if the flood impacts the flow resistance. Thus, to account for the floodplain drag when h>h_f, the Manning coefficient is reevaluated as follows:

$$\begin{array}{}\text{(28)}& {n_{{\mathrm{cf}}_i}}^t={\displaystyle \frac{n_i}{{{\mathit{\zeta }}_{n_i}}^t}}.\end{array}$$

In the calibrated model, C_nfp values along the mainstem range between 0.3 and 1. Under typical flood conditions, Y varies between 0 and about 0.3, resulting in an increase of the effective Manning coefficient of approximately 0 %–6 % (see S5 in the Supplement). Larger corrections, up to about 20 %, may occur only during extreme simulated floods when Y approaches 0.5. Such variations are consistent with resistance changes reported in studies of compound channel hydraulics and floodplain–channel momentum exchange (e.g., Nicollet and Uan, 1979; Smart, 1992; Knight and Shiono, 1996; Bousmar and Zech, 1999). The coefficient C_nfp was calibrated during the model calibration phase to reproduce observed stage–discharge relationships at the mainstem gauging stations.

Following the same reasoning as for velocities the ratio of the dimensionless grain-related bed shear stresses ${\mathit{\theta }}_{\mathrm{c}}^{\prime }$ and ${\mathit{\theta }}_{\mathrm{cf}}^{\prime }$ should also be a function of ζ_n. Indeed, according to Eq. (14):

$$\begin{array}{}\text{(29)}& {\displaystyle \frac{{\mathit{\theta }}_{\mathrm{cf}}^{\prime }}{{\mathit{\theta }}_{\mathrm{c}}^{\prime }}}={\left({\displaystyle \frac{f_{\mathrm{cf}}^{\prime }{u}_{\mathrm{cf}}}{f_{\mathrm{c}}^{\prime }{u}_{\mathrm{c}}}}\right)}^{\mathrm{2}}={\mathit{\zeta }}_n^{\mathrm{2}}.\end{array}$$

Here, $f_{\mathrm{cf}}^{\prime }=f_{\mathrm{c}}^{\prime }$ is assumed, as f^′ is a grain-related friction factor, not a flow resistance factor (Yen, 2002): the floodplain drag is already accounted for in u_cf, through n_cf. The shear velocity term used for calculating the Rouse number (Eq. 10) should be also affected by the floodplain drag:

$$\begin{array}{}\text{(30)}& u_{\ast \mathrm{cf}}={\mathit{\zeta }}_n{u}_{\ast \mathrm{c}}.\end{array}$$

However, when shifting from a 1D to a 2D framework, the transverse profiles of ${\mathit{\theta }}^{{}^{\prime }}$, and consequently of $C_{\mathrm{b}}^{\ast }$ and $Q_{{\mathrm{s}}_{\mathrm{cf}}}^{\ast }$, are likely to be more strongly affected than the lateral profile of the depth-averaged velocity (see S6 in the Supplement), in particular near the banks. To account for the complex 2D effects on sediment transport capacity, effects not considered in the initial computation of the transport capacity $Q_{{\mathrm{s}}_{\mathrm{cf}\mathrm{0}}}^{\ast }$ which was initially calculated using the corrections for ${\mathit{\theta }}^{{}^{\prime }}$ and u_∗ corrections in Eqs. (29) and (30), the following formulation is applied when the floodplain is active (i.e. when h>h_f):

$$\begin{array}{}\text{(31)}& {{Q_{\mathrm{s}}^{\ast }}_{{\mathrm{cf}}_i}}^t={{\mathit{\zeta }}_n^{\mathit{\eta }}}_i^t{{Q_{\mathrm{s}}^{\ast }}_{\mathrm{cf}{\mathrm{0}}_i}}^t,\end{array}$$

with η an exponent to calibrate which accounts for these complex 2D effects.

4.3.2 Bed roughness influence for low waters

In the large Amazonian rivers, a decrease in bed roughness influence with increasing water levels has been observed (see example in S7 in the Supplement). To model this in the SWAT-Amazon version, the Manning coefficient is modified using the factor ζ_nch, defined as:

$$\begin{array}{}\text{(32)}& {{\mathit{\zeta }}_{n{\mathrm{ch}}_i}}^t=\mathrm{1}+C_{n{\mathrm{ch}}_i}\left({\displaystyle \frac{h_{{\mathrm{ch}}_i}-h_i^t}{h_i^t}}\right),\end{array}$$

where C_nch (–) is a coefficient and h_ch (m) the water height below which the additional bed roughness correction activates. Calibrated values of h_ch range between 35 % and 75 % of the maximum simulated water level across the mainstem reaches, and C_nch between 0.12 and 0.4, resulting in an increase of the effective Manning coefficient of approximately 10 %–30 % at low stage, decreasing progressively as stage rises above h_ch (Fig. S4). Both parameters were calibrated to reproduce the observed stage–discharge and stage–velocity relationships at the mainstem gauging stations.

4.4 Calibration and sensitivity analysis

The model calibration was performed using the SWATrunR package (Schürz, 2019), which enables parallel processing. To run the SWAT-Amazon executable and calibrate parameters, including the newly introduced ones (Table 1), an R-Notebook was written (Fig. 1d). It allows users to export interactive figures and perform sensitivity analyses. Both SWAT-Amazon and its R-Notebook for calibration are available for download at: https://github.com/william-santini/SWAT-Amazon (last access: 6 June 2025).

5 Results

5.1 Water discharge simulations

For water discharge, model calibration was performed over the 2010–2015 period based on multiple hydraulic diagnostics (water levels, velocities, and stage–discharge relationships). Model performance was further evaluated using independent hold-out periods and direct comparisons with gauging measurements that bypass rating curves, as detailed in S9 in the Supplement. These complementary evaluations show consistent performance across periods and support the temporal robustness of the model.

At a daily time step, SWAT-Amazon simulations at the watershed outlet show excellent performance (Fig. 3d): NSE (Nash-Sutcliffe Efficiency) = 0.92, KGE (Kling–Gupta efficiency) = 0.95, PBIAS (Percent Bias) = −1.8 %, LogNSE (NSE on the logarithms of the series) = 0.92 over the 2000–2016 observation period used here for inter-station evaluation (see Moriasi et al. (2007) for details on these metrics). SWAT-Amazon significantly improves over the standard SWAT model (NSE = 0.86), using Muskingum routing with maximum flood attenuation. Moreover, the standard SWAT simulation predicts flood peaks 1–2 months earlier than observed, whereas SWAT-Amazon correctly synchronizes them. The model accurately captures hydrological dynamics and interannual variability. As highlighted by Yamazaki et al. (2011), the difference between kinematic and diffusive wave simulations was minimal (not shown), confirming that Ucayali flood attenuation mainly results from floodplain buffering.

Figure 3 Water discharge simulations for the stations with gauging data. (a–d) Observed discharge (marine blue) vs. simulated discharge using the default SWAT model (orange) and SWAT-Amazon (green) at a daily time step, with punctual ADCP gauging values (blue circles). (e) Full 37-year simulation (1983–2019) with SWAT-Amazon for the same station group.

At Pucallpa (Fig. 3c), SWAT-Amazon shows only slight improvements over standard SWAT due to the less developed floodplain. At the Andean outlet (Lagarto and Puerto Inca), where floodplain influence is minimal, both models perform similarly, though SWAT-Amazon slightly outperforms the default version. Despite a good daily NSE (0.72) at Puerto Inca, the model struggles to reproduce rapid flood oscillations typical of piedmont hydrographs. This issue, independent of the routing model, stems from uncertainties in rainfall estimation. Before final calibration, systematic biases (−20 % to +20 %) were observed, with underestimation in piedmont stations and overestimation in plains, primarily due to the precipitation dataset. These biases were corrected using interannual adjustment factors in SWAT .sub files. Additional errors in the PISCO dataset were identified and corrected by standardizing precipitation time series across station subgroups. However, PISCO still underestimated precipitation between Contamana and Requena for 2016, 2017, and 2019, leading to their exclusion from efficiency calculations.

Despite these limitations, the bias-corrected PISCO dataset demonstrated a high degree of homogeneity and robustness, allowing extension of observations across all stations (virtual and conventional) for simulations covering 1983–2019 (Fig. 4e), adding 13 years at Requena and 26 years at Lagarto. This extension is particularly valuable for future studies in this poorly monitored region, especially given the high accuracy of the weighted seasonal Nash-Sutcliffe Efficiency for low (wsNSE_LF) and high (wsNSE_HF) flows (see Zambrano-Bigiarini and Bellin, 2012, for details).

Figure 4 Example of (a) water levels and (b) velocities simulation at Pucallpa. Marine blue line: observations, orange line: best default SWAT simulation with Muskingum, green line: SWAT-Amazon simulation, blue filled circle: ADCP gauging values. Error bars, which were less than 3 % for ADCP measurements, are not shown for clarity.

5.2 Water levels, velocities, and rating curves simulations

The hydraulic radius R_h was used to compare simulations with observations (Figs. 4 and 5). Indeed, for large cross-sections, the mean depth h_m (m) approximates R_h, and for the modelled rectangular cross-sections, $R_{\mathrm{h}}\approx h=h_{\mathrm{m}}$. Moreover, observed water levels were not directly comparable due to offset differences in staff gauge zero-values relative to assumed river bottom elevation. In standard SWAT, once the bankfull height is exceeded, flow instantly spreads into the floodplain, forming a single cross-section instead of the usual approach in hydraulics of separating channel and floodplain flows (e.g. Einstein, 1950; Einstein and Barbarossa, 1952; Yen, 2002). This sudden change in cross-section geometry thus causes a discontinuity in flow velocity and hydraulic radius (Fig. 4), since P_h increases sharply while A_h grows more moderately ($R_{\mathrm{h}}=A_{\mathrm{h}}/P_{\mathrm{h}}$). Below bankfull height, standard Muskingum simulations overestimate hydraulic radius and underestimate velocities (Fig. 4), and it is not possible to calibrate u(h) rating curves. Therefore, the standard SWAT model is not able to simulate realistic hydraulic radius (water levels) and velocities and even less sand loads with transport capacity laws, for which these variables are required. Conversely, SWAT-Amazon generates robust daily water level and velocity time series, closely matching observations (Fig. 4), with NSE values between 0.77 and 0.93 for water levels and 0.79 to 0.92 for velocities, the lowest at Puerto Inca, while all others exceed 0.89. It produces consistent Q(R_h) and u(R_h) rating curves (Fig. 5), accurately capturing slope-controlled hysteresis and “duckbill” damping when h_f is exceeded, as Manning's coefficient increases with relative water height due to floodplain effects.

Figure 5 Rating curves between hydraulic radius (R_h) and discharge (left) or velocity (right). Green circles: SWAT-Amazon simulation, blue filled circles: ADCP measurements.

5.3 Sand routing

Calibration focused on the September 2009–Augsut 2015 period, when sediment monitoring protocols were enhanced, including higher sampling frequency at Requena between November 2012 and June 2013, where one sample was collected each 2 d plus three sampling repetitions each 10 d. Beyond, sampling was conducted at 5 d intervals during the wet period between July 2013 and September 2015. Additionally, the concentration gaugings were performed in all sites with a higher number of samples collected throughout the cross-section, particularly in the first half of the water column, to ensure more accurate sand concentration calculations. Outside this interval, uncertainties in sand flux observations increase, which complicates the definition of a robust and fully independent validation period (see S9 in the Supplement for a summary of model performance across periods). The lower performance outside the calibration period primarily reflects uncertainties in rainfall forcing and observations, particularly during rapid Andean flood events, rather than a degradation of model performance. To further evaluate model performance, direct comparisons between simulations and gauging measurements were performed (Fig. S6). These complementary evaluations support the temporal robustness of the model despite observational limitations and data heterogeneity.

At Lagarto, the sand routing model accurately reproduces sand fluxes (Fig. 6a, daily NSE = 0.8), validating the capacity-limited flux assumption at the Andean outlet (cf. Sect. 3.4). At Puerto Inca (Fig. 6b, NSE = 0.44), the model struggles due to rainfall data and sharp flux peaks. Nevertheless, at the mainstem level of the Ucayali, the influence of this discrepancy is limited. At Requena, the model closely matches observed sand fluxes (Fig. 6c, NSE = 0.86) with peaks coinciding with maximum rainfall in January-February. From March, sand flux decreases while flow increases, indicating no correlation between sand flux and discharge. This decline is concomitant with the crossing of the threshold h_f from which the floodplain watering impacts the transport capacity. In the 2010 drought year, the river briefly reached this threshold, with minimal impact on sand flux: Q and Q_ss are well correlated. In the remaining years, a second sand peak occurred in May–June. This is concurrent with the recovery of river transport capacity, which is enabled by a reduction in flow resistance due to the dewatering of the floodplain and an increase in energy availability in the main channel, resulting from the influx of floodplain and black water supplies, which have low sediment concentrations. The 2012 extreme flood event, intensively monitored, highlights this key process for sediment routing dynamics.

Figure 6 Observed and simulated sand fluxes for the gauging stations with sediment monitoring: (a) Lagarto, (b) Puerto Inca, (c) Requena. Gray line with stars: observations, green line: SWAT-Amazon simulations, orange filled circles with error bars: gauged sand flux values. Observed (blue dashed line) and simulated (green dashed line) water discharge at Requena are plotted in (c). The blue horizontal line represents the discharge triggering floodplain activation (corresponding to h_f), approximated due to the non-bijective stage-discharge relationship.

5.4 Remote-sensed fine sediments fluxes

Calibration campaigns established a relationship between fine sediment concentration C_f (mg L⁻¹) (at z=h) and the NIR-to-red reflectance ratio R, supporting a single model for the entire basin (Fig. 7) with a high coefficient of determination (R² = 0.94) and a low Mean Absolute Error (MAE = 59 mg L⁻¹):

$$\begin{array}{}\text{(33)}& C_{\mathrm{f}}\left(h\right)=\mathrm{822.8}\times R^{\mathrm{2.86}}+\mathrm{0.001}\times \exp \left(\mathrm{10.4}\times R\right).\end{array}$$

This model accounts for reflectance saturation in the red band at high concentrations, providing a better fit across the full concentration range than a simple power-law equation. It was validated across all hydrological conditions from 2000 to 2019 using matchups between time series of in situ fine sediment concentrations monitored at Requena and Lagarto and co-located satellite reflectance ratios at a monthly time step (R² = 0.78, MAE = 132 mg L⁻¹) (Fig. 7). Note that Eq. (33) is already corrected for adjacency effects, through a simple offset of +0.2 applied to the reflectance ratio to account for water pixel contamination by riverbanks. Finally, the calibration campaigns confirmed that the contribution of suspended sand to the reflectance ratio R is negligible (see S8 in the Supplement).

Figure 7 Relationship between fine sediment concentration (C_f) at the water surface and the ratio R of NIR to red reflectance. Orange dots represent calibration points, based on 42 field measurements in the Ucayali Basin, where reflectance was measured using a hyperspectral radiometer and fine sediment concentrations were measured at the water surface. Gray dots correspond to matchups between satellite-derived reflectance ratio and fine sediment concentrations monitored at gauging stations, averaged at a monthly time step.

5.5 Validation of the integrated approach at super stations

The integration of remote sensing and hydrological modelling was validated at two super stations in the basin: Requena (Fig. 8) and Lagarto (not shown). Total sediment fluxes (Fig. 8c) were calculated by summing: (i) sand flux (Fig. 8a) from SWAT-Amazon simulations and (ii) fine sediment flux (Fig. 8b) from satellite data combined with SWAT-Amazon flows. The results align well with in situ flux measurements (daily NSE: 0.87 at Requena, 0.79 at Lagarto, monthly NSE: 0.87 at Requena, 0.86 at Lagarto), and suggests that both stations could be monitored in this way with a few calibration campaigns. This supports the definition of minimal observational requirements for transferring the method to other Amazonian basins (see Sect. 6.2.3). The evaluation focuses on variables that directly constrain main-channel hydrology and sediment transport, and thereby indirectly constrain floodplain processes. In the present framework, floodplains are represented using a simplified storage approach, in which water levels are spatially aggregated at the reach scale. However, natural floodplains exhibit strong spatial heterogeneity, including topographic gradients, secondary channels, lakes, and varying degrees of hydraulic connectivity. As a result, local water levels may not directly reflect the large-scale storage dynamics represented in the model. In addition, existing satellite products lack the spatial and temporal resolution required to robustly constrain floodplain water levels and associated storage volumes, particularly given the diversity of water surface types (e.g. secondary channels, connected lakes and isolated water bodies) and their distinct behaviours. Accordingly, local floodplain water levels are not explicitly reproduced; instead, the analysis focuses on constraining water and sediment fluxes and associated storage volumes at the subbasin scale using multi-source observations of main-channel dynamics. The increasing availability of satellite observations over time (Fig. 8b) further strengthens the observational constraints on these fluxes, with missing data mainly occurring at the beginning of the time series due to the progressive deployment of satellite sensors (Terra since 1999, Aqua since 2002, VIIRS since 2012).

Figure 8 Integrated monitoring of sediment fluxes at the basin outlet. The gray line with stars shows observed sediment fluxes, the green line represents SWAT-Amazon simulations, and orange filled circles indicate gauged flux values.

6 Discussion

This study presents the first integrated approach for monitoring hydro-sedimentary fluxes in the Amazon basin, providing the most extensive and robust daily time series of hydro-sedimentary fluxes for the Upper Amazon (20 years for fine sediments to 37 years for water discharge and sand flux) with exceptionally high NSE values. It introduces a physically based methodology for modelling transport capacity, using realistic hydraulic parameters derived from the calibration of u(h) and Q(h) rating curves.

Previous large-scale efforts with the MGB (Modelo de Grandes Bacias) model (Fagundes et al., 2021, 2023), contributed significantly to understanding sediment transport across South America. However, challenges remain in representing sand transport, particularly its suspension dynamics. The MGB model assumes that sand transport is predominantly bedload, whereas field observations indicate that sand can account for a substantial fraction of the suspended sediment load in large Amazonian rivers, reaching up to ∼ 70 %. Rouse numbers in the range 0.2–0.8 further indicate transport in graded suspension rather than intermittent transport (Santini et al., 2019; Martinelli, 2022). This conceptual difference likely leads to underestimates of sediment load in the Ucayali Basin in previous MGB applications, with reported estimates being up to nearly three times lower than observations, although its contribution would need to be confirmed through a controlled inter-model comparison. These findings underscore the necessity of detailed, regionally focused studies based on long-term measurements and rigorous data consistency analyses, rather than broad continental-scale assessments, which often involve considerable uncertainties and may lead to misinformed sustainable development strategies and mitigation policies.

For the first time, satellite-based sediment monitoring is applied exclusively to the fine fraction. A relationship (Eq. 33) is proposed that is independent of hydrodynamic fluctuations affecting the surface concentration, since sand contributes from ∼ 5 % to 50 % to the surface concentration in the Ucayali. This contrasts with previous studies (Espinoza-Villar et al., 2012, 2013, 2017; Park and Latrubesse, 2014; Martinez et al., 2015) where satellite reflectance was solely compared with total sediment concentration or where remote-sensing was only used to calibrate the model (Fagundes et al., 2020).

However, Pinet (2017) noted hysteresis in Madeira River relationships due to variations in grain diameter at the water surface, and Santini (2020) suggested that a fraction of fine sediments might also be sensitive to turbulence-induced lift variations. This may introduce limited variability in the relationship during specific hydrological phases, particularly resuspension events and low-water conditions. However, previous studies have shown that such effects are strongly reduced when using the spectral band ratio (NIR $/$ RED) (Martinez et al., 2015; Pinet, 2017), with remaining variability comparable to satellite reflectance uncertainty. In addition, evaluation against satellite–in situ matchups across a wide range of hydrological conditions does not reveal any systematic hysteresis at the monthly time scale. As a result, the impact on long-term sediment flux estimates is expected to remain limited relative to other sources of uncertainty and to the overall magnitude of sediment transport.

Before drawing conclusions and interpreting the mass balances (Sect. 6.4), it is crucial to assess the robustness and limitations of the method to ensure that the necessary nuances are applied, particularly at the virtual stations. To this end, a sensitivity analysis of the SWAT-Amazon model was conducted.

6.1 Model sensitivity and equifinality

The Sobol variance-based sensitivity analysis serves a diagnostic purpose: it identifies which parameters require precise measurements for SWAT-Amazon application and which can be assigned wider priors without degrading budget estimates. This hierarchy directly informs the implementation of the Generalized Likelihood Uncertainty Estimation (GLUE) framework (Beven and Binley, 1992), used here for uncertainty propagation (Sect. 6.3.2 and S12 in the Supplement).

The Sobol method was applied to sub-basin 21 over the 2009–2015 period to assess the relative importance of SWAT-Amazon input parameters and their interactions, identifying uncertainties from the model's formalisms. First-order indices (S_i) and total-order indices (T_i) were computed to quantify each parameter's individual contribution (with S_i) and its total effect (with T_i), including interactions with other parameters, on the overall variance of the model output. The NSE criterion was applied for global outputs (Fig. 9a, d), and the interannual daily average for time-series outputs (Fig. 9c, f). The analysis focused on a selected set of SWAT-Amazon parameters, based on calibration experience and physical understanding, rather than exhaustively testing all parameters. Wide prior ranges were deliberately used to ensure a comprehensive exploration of the parameter space (S10 in the Supplement).

Figure 9 Sobol sensitivity analysis results for water (left) and sand routing (right), applied to sub-basin 21 over the 2009–2015 period. (a, d) The analysis based on the NSE criterion, (c, f) summary of the temporal analysis as interannual daily averages. (b, e) The interannual daily average of the observed flux (blue line), the SWAT-Amazon simulation (green line), and the envelope encompassing all simulations performed for the Sobol analysis (light-green ribbon). The blue horizontal line on (b) represents the approximative discharge triggering floodplain activation.

6.1.1 Water routing

The analysis focused on the parameters set (n, B, h_f, C_nfp, k_fp), with n accounting for all the resistance term $S_{\mathrm{b}}/n$. Results show a greater sensitivity of the model to (h_f, C_nfp, k_fp), which control flood wave attenuation by the floodplain (Fig. 9a). The interannual time-series analysis (Fig. 9c) reveals greater sensitivity during flood recession than rising waters. The recession is particularly challenging to calibrate due to rapid water discharge drops, where even slight shifts in the timing of floodplain recession cause large discrepancies between simulated and observed flows throughout the entire recession limb. Therefore, the parameter h_f, controlling recession onset, must be carefully assessed. The strong oscillation of Sobol indices for n during low waters (Fig. 9c) reflects the impact of h (through B) on bed roughness and flow resistance when h<h_ch (cf. Eq. 32). Small variations of h induce large changes in n during low waters, but with minimal impact on discharge, as shown by the interannual discharge plot in Fig. 9b. This effect diminishes quickly as water levels rise.

6.1.2 Sand routing

The sensitivity analysis, using the parameter set (β_s, C_bk, d_b, d_s, η, k_bed), shows that d_b is the most critical calibration parameter (Fig. 9d, f), as previously highlighted by Fagundes et al. (2023). In the present framework, d_b and d_s are assumed to remain constant over the seasonal cycle. Calibrated d_b values for sub-basins 19, 5, 14, 23, 22, and 21 are 252, 240, 240, 242, 235, and 220 µm, respectively, matching observed mean diameters at Lagarto (260 µm), Pucallpa (243 µm), and Requena (228 µm). Calibrated d_s values are approximately 80 µm for all sub-basins, except Lagarto (98 µm), are consistent with PSD observations. Since the calibration of the flood recession limb directly influences sand resuspension (Eq. 22), h_f emerges also as a key parameter in sand routing. The influence of k_bed is minimal, as the sand load input in sub-basin 21 is sufficient, eliminating the need for riverbed erosion to compensate for sediment supply deficits.

6.1.3 Equifinality analysis

GLUE-based dotty plots were produced by sampling the Sobol-influential parameters over wide ranges (S11 in the Supplement). The water routing parameters (h_f, B, k_fp, C_nfp, n) all show identifiable optima to varying degrees depending on the reach, consistent with the Sobol analysis. Critically, d_b exhibits a physical optimum near 220–250 µm matching observed grain sizes, and spurious coarser modes, precluding sand suspension. This bimodality highlights a well-known limitation of purely statistical calibration: without physical constraints, optimization can converge toward parameter values that are mathematically acceptable but physically meaningless. Overall, since d_b is well constrained and water routing parameters are collectively identifiable, the hydro-sedimentary budgets remain robustly constrained despite residual equifinality across all behavioural simulations (S12 in the Supplement). These findings justify the use of expert measurement-informed priors for uncertainty propagation.

6.2 Calibration insights from the SWAT-Amazon method

In the hydraulic routing method (Sect. 4.1), compensating hydraulic effects limit the influence of B and n on routed discharge. Reducing n increases flow velocity u and discharge Q (Eq. 5), which decreases the water volume stored in the reach and therefore h (Eq. 1). The resulting decrease in cross-sectional area leads to a proportional reduction in u and Q.

Therefore, when h<h_f, discharge calibration relies exclusively on the default SWAT rainfall-runoff model. However, u and h are directly related to the parametrization of n and B and are the primary variables in the transport capacity model utilized here (Sect. 4.2). Q is also a key variable (Eq. 9). As B is poorly-known, it was excluded from Eqs. (6) and (8) of the transport capacity model by replacing B h with $Q/u$ and u B with $Q/h$. With regard to n, another significant source of uncertainty, it should be noted that this variable is not included in the transport capacity equations. This is because u is used and calibrated prior to the Q_s computation.

6.2.1 Calibration strategy for a super station

According to the previous analyses, the calibration strategy for stations with robust, long-term hydro-sedimentary monitoring is:

• a.

  Start by calibrating Q in each reach, for h_f only, using the SWAT's default hydrologic parameters.

• b.

  Calibrate Q, considering floodplain effects, using h_f, C_nfp and k_fp (or θ_fp).

• c.

  Calibrate u and h by adjusting n and B only; Q is unaffected by this calibration.

• d.

  Check the relationships Q(h) and u(h), revisiting step c if needed.

• e.

  Compute the $Q_{\mathrm{s}}(h,u,Q)$, independently of n and B. If necessary, adjust Q_s using parameters in Table 1, particularly d_b, the most sensitive parameter.

It is important to emphasize that the optimal calibration for water discharge may not align with the best calibration for water level, velocity, and sand load time series. A compromise must be made.

6.2.2 Calibrating virtual and low-data stations with satellite altimetry

The sensitivity of the model to n and B emphasizes the importance of h (Eqs. 1, 2, 4) in the routing module, particularly for floodplain inundation and dewatering triggers (Sect. 4.1), as already discussed above. For ungauged reaches (5, 23, 22, see Fig. 2a) in the plain, n and B can be calibrated using relative water levels (hydraulic radius) derived from spatial altimetry and shifted by a constant offset to match simulations (Fig. 10). The Manning coefficient is supposed already known to a reasonable degree, as coherent values of $\mathrm{1}/\mathrm{45}$ ± 0.0005 s m${}^{-\mathrm{1}/\mathrm{3}}$ have been found across all the other plain's sub-basins. Given the large sub-basins and significant channel lengths, bed slope errors are assumed to be smoothed, and therefore S_b also well assessed. Thus, velocities can be reliably determined if B is calibrated through h calibration. The calibration for reaches (5, 23, 22) was also guided by looking for the best simulation at Pucallpa and Requena.

Figure 10 Calibration of ungauged reaches using satellite altimetry: Satellite-derived water levels (blue points) shifted by a constant offset to match simulations: triangles (Envisat), circles (Jason-2), squares (Jason-3), diamonds (Sentinel-3B). Conventional in situ observations appear as an orange line, while the curated series incorporating satellite altimetry is in blue. SWAT-Amazon simulations are shown in green. (a) Reach 5 (virtual station S3B331-R5), (b) Reach 23 (Contamana), (c) Reach 22 (virtual station JA204-S3B310-R22).

Finally, implementing an integrated approach combining modelling and remote sensing for monitoring hydro-sedimentary fluxes would require revisiting hydrological network operations. Specifically, a few gauging campaigns (around four, strategically timed along the annual hydrograph) could be sufficient for model calibration. Additionally, synchronizing gauging with satellite altimetry would enable the construction and simulation of robust rating curves.

6.2.3 Minimum observational requirements

Based on the long-term monitoring experience of the CZO HyBAm network and the results obtained in this study, the framework can be applied to other Amazonian basins provided that a minimal observational dataset is available at a limited number of super stations and complementary virtual stations. For the last ones, this includes: (i) approximately four to eight discharge gauging campaigns covering the main stages of the annual hydrograph (rising limb, high water, falling limb, and low water), (ii) a comparable number of suspended sediment gaugings distributed along the hydrograph and carried out using standardized and hydraulic-based protocols (Santini et al., 2019), (iii) one to two detailed PSD surveys in relatively homogeneous upstream reaches, increasing to four to five in more heterogeneous downstream settings (e.g. below major confluences), and (iv) two to three radiometric calibration campaigns spanning contrasted hydrological conditions.

6.3 Framework limitations and uncertainty

6.3.1 Model structural limitations

This study adopts a hydraulic representation of the river, modelling the floodplain as a reservoir. Rapid exchanges through secondary channels are neglected, and their influence is approximated on discharge propagation through a minor adjustment of the Manning's coefficient. Guilhen et al. (2022) used a previous version of SWAT-Amazon with a slightly modified floodplain configuration, allowing 1d floodplain flow through its whole cross-section with a Manning-Strickler equation. However, in addition to being almost impossible to calibrate, the floodplain flow being extremely sensitive to the floodplain Manning coefficient because of its huge cross-section, the changes in terms of results were extremely negligible: less than 1 % of the water discharge in the main channel, well below expected model precision.

Another key assumption in SWAT-Amazon is the absence of backflow from the floodplain to the main channel during the rising limb. Observations in the Amazon basin support this, as low-gradient floodplain channels often experience tributary flow blockage, with the Tapajós River in Brazil being an extreme example. However, during deflooding, floodplain drainage dynamics differ from flooding processes. In SWAT-Amazon, these dynamics are governed by floodplain cross-section geometry. Introducing a Maillet-type law could improve recession representation, but discrepancies in timing appear mainly driven by rainfall variability, masking this effect. Thus, no modifications were made.

Floodplain sediment concentration during deflooding was assumed constant. Nevertheless, the unique monitoring of sediment concentration in a floodplain channel, performed at Lago Grande de Curuai (Brazilian Amazon) shows a decrease from ∼ 800 mg L⁻¹ during low-waters when the floodplain water flows back into the main channel to ∼ 40 mg L⁻¹ when the channel is under the influence of the Amazon River, with a mean concentration of 135 mg L⁻¹ (Moreira-Turcq et al., 2013). This value appears to be relatively close to that calibrated for C_bk (for the sand fraction only) at sub-basins 23, 22, and 21 (∼ 180 mg L⁻¹).

A dynamic law linking C_bk to water level could refine peak resuspension modelling but requires concentration data for calibration. Alternatively, Fagundes et al. (2023) treated the floodplain as a homogeneous sediment reservoir, balancing suspended fluxes with settling. However, this method applies only to clay and silt, excluding sand, and the settling law can be challenging to calibrate in the absence of observational data. Given floodplain water storage timescales (months), most fine sediments likely settle before being remobilized. Furthermore, as underscored in the introduction, sediments resuspended during recession can be centuries to millennia old, indicating long-term floodplain recycling rather than short-term remobilization.

Lastly, the distinction between E_bk and E_bed is partly supported by the calibration experiments. When bed erosion is activated (k_bed>0), the model tends to generate rapid and abrupt peaks in simulated sand flux once the transport capacity exceeds the available sand load, whereas adjusting C_bk produces a smoother and more progressive increase in sand concentration that better reproduces the secondary peaks observed during the recession phase. The calibrated values of K_bed remain very small (K_bed≪1), suggesting that riverbed erosion contributes only weakly to the simulated suspended sand flux in the main stem, although a minor contribution cannot be fully excluded.

6.3.2 Framework uncertainty

A GLUE-based uncertainty analysis was conducted over a Latin Hypercube ensemble of 2500 simulations (S12 in the Supplement), scoped to routing and sand transport parameters, conditional on the hydrological forcing, with physically informed priors (Sect. 6.1).

For discharge and sand flux, GLUE-weighted envelopes (5th–95th percentiles) remain narrow relative to the interannual amplitude, indicating that parameter uncertainty propagates weakly compared to the magnitude of hydrological variability and that the dominant dynamics are robustly constrained. Envelope width depends jointly on data availability and measurement quality, which act as coupled controls on parameter identifiability rather than independent sources of uncertainty. Floodplain water storage uncertainty in particular is strongly controlled by h_f, emphasizing the need for accurate stage–discharge rating curves to constrain this parameter reliably (Sect. 6.2.3). Importantly, because the budget indicators (trapping and recycling fractions) depend primarily on h_f, B, and d_b, which are well-constrained parameters (Sect. 6.1.3), they remain robust across the behavioural ensemble despite the wider prior ranges assigned to less identifiable floodplain parameters. This structural stability is confirmed by the GLUE-propagated budget envelopes (S11 and S12 in the Supplement).

Fine sediment fluxes achieve low bias and good performance metrics (Fig. 8b), yet exhibit wider uncertainty envelopes. This apparent performance is largely driven by the dominant seasonal signal rather than true predictive skill. Satellite retrieval errors cancel over long-term integration, masking short-term variability: the method reliably captures the mean seasonal cycle, but not its anomalies. Satellite radiometry is therefore robust at the climatological scale but not at the event scale.

Both conclusions rest on a common foundation: physically meaningful uncertainty quantification requires prior knowledge, which itself requires long-term, multi-variable observations, as provided by networks such as the CZO HyBAm. Uncertainty is therefore not only a modelling issue, but fundamentally a data-structure constraint. Model calibration and interpretation depend on sustained fieldwork and on a detailed understanding of measurement protocols and data limitations. As a result, rigorous data governance, FAIR practices, and strong research–operational links are essential for reproducibility and long-term continuity.

6.4 Hydro-sedimentary dynamics in Amazonian foreland systems: long-term insights

The robustness of the results enables precise water and sediment balances, distinguishing fine particles from sand, at each virtual and conventional station (Fig. 11). This is achieved at an unprecedented spatiotemporal resolution over the periods 1983–2019 (37 years) for water and sand fluxes and 2000–2019 (20 years) for fine sediments.

Figure 11 Interannual water (1983–2019) and sediment (2000–2019) balances along the Ucayali River. Lines show interannual monthly means and shaded bands represent the 5th–95th percentile envelope of the behavioural GLUE ensemble (2500 runs). (a) Water Discharge. (b) Floodplain water storage (V_fp) in km³. (c) Normalization of V_fp by the reach's length (Δx) for cross-sub-basin comparison. (d) Fines suspended sediments. (e) Suspended sand fraction. (f) Total suspended sediment load. See S12 in the Supplement for full uncertainty propagation details.

6.4.1 Dynamics of flood waves, flooding, and sediment transport

Water and sediment peak fluxes at the Andean outlet occur in February (Fig. 11), declining sharply from March to May as Andean precipitation decreases. The sediment flux at the basin outlet is strongly correlated with Andean discharge (Fig. 12c), confirming the Andes as the primary sediment source. This results from intense erosion along the Eastern Cordillera and Sub-Andean thrust belt, where high precipitation erodes Paleozoic and Mesozoic sedimentary sequences incorporated into Cenozoic deposits eastward as the orogenic prism advances through crustal shortening and thickening (McQuarrie, 2002a, b; Espurt et al., 2008; Gautheron et al., 2013; Pfiffner and Gonzalez, 2013). Additionally, the Central Andes' convex hypsometric profile (Montgomery et al., 2001) promotes deep fluvial incision.

Figure 12 Identification of sediment sources at the interannual scale. (a) Water discharge components: Andean Basin outflow (Q_Andes), basin outlet discharge excluding lateral contributions (ΔQ_plain) and floodplain flows (Q_Floodplain), all in m³ s⁻¹. (b) Decomposition of total sediment flux at the basin outlet Q_t,21 (10⁶ t d⁻¹) by sediment sources: red for Andean erosion, blue for floodplain recycling. (c) Relationship between Q_t,21 and ΔQ_plain (circles) and Q_Andes (diamonds). The red regression line represents the relationship between Q_Andes and Q_t,21 during rising limb months (d) Relationship between Q_t,21 and Q_Floodplain. The blue regression line illustrates the relationship between Q_Andes and Q_t,21 during recession months.

At the lowland outlet, discharge peaks two months after the Andes (April–May) delayed primarily by alluvial plain storage, diffusive flood wave propagation, and runoff from the northern part of the basin, where seasonality, driven by the South American monsoon system (Garreaud et al., 2009), is less pronounced than in the south. This upstream-to-downstream flooding dynamic leads to progressive floodplain inundation (February–April, Fig. 11b, c) and gradual drainage (March–July). Floodplain backflow (Fig. 12a) drives sediment remobilization from April to August (Fig. 12b, d). This secondary sediment source accounts for ∼ 22 % of the total suspended sediment load at the basin outlet and can have significant impacts on river dynamics, with pronounced effects during major floods (e.g., 2012, Fig. 6c). Floodplain flows, comprising precipitation, river inflows and backflows, are more significant than the water entering the floodplain from the main river (Fig. 12a), reinforcing the assumption that floodplain waters are blocked by the main channel during the rising limb.

In the following analysis, the fine sediment load of the Pachitea Basin was estimated using a rating curve derived from the relationship between water discharge and measured fine sediment load, based on simulated flows. For the other Andean sub-basins (8, 9, 3), fine sediment load was regionalized according to the drainage area of their respective Andean catchments. The Ucayali floodplain consists of three geomorphological compartments (Fig. 2b), each characterized by distinct processes:

Lagarto – Pucallpa (Box 1, Fig. 2b): A major thrust fault system involving the crystalline basement exerts significant tectonic control over the narrow sedimentary basin wedged between the young Shira Mountains and the Fitzcarrald Arch (Espurt et al., 2008; Pfiffner and Gonzalez, 2013; Gautheron et al., 2013). Here, 36 % of the 420 [386, 454] × 10⁶ t yr⁻¹ of Andean sediment accumulates, with sand deposition (−100 [−83, −117] × 10⁶ t yr⁻¹ ) accounting for 66 %. This results in alluvial plain aggradation at 3.7 mm yr⁻¹ (9750 t km⁻² yr⁻¹), highlighting strong neotectonic control. In this compartment, water storage reaches 1.5 [1.1, 1.9] km³ in March on average, with flood durations (days with h ≥ h_f) varying from 0 to 56 d in sub-basins 19 and 5, respectively. Consequently, daily-scale Andean flash floods superimposed on the annual monomodal flood wave are rapidly attenuated, while downstream, water storage increases as the floodplain expands. This process also dampens high-amplitude fine sediment peaks (−51 [−66, −35] × 10⁶ t yr⁻¹) through the activation of flood channels. These channels, remnants of the river's past mobility, further enhance the buffering effect. The sudden floodplain expansion at sub-basin 14 as the river course is released from the tectonic constraint imposed by the Shira Mountains reduces fine sediment loads (Fig. 11d), while high meander migration rates (Li et al., 2023) indicate active lateral erosion and deposition.

Pucallpa – Contaya Arch (Box 2, Fig. 2b). This section is a large piggyback basin located between the Andean wedge top and the Moa Divisor thrust fault system, with increased floodplain storage and longer flood durations (78–125 d). In March, it retains 10.4 [8.3, 12.5] km³ of water (∼ 20 % of discharge). Fine sediment load increases by +10 [+2, +17] × 10⁶ t yr⁻¹ due to Andean inputs, mainly from the Pachitea sub-basin, where one of the major Andean precipitation hotspots drives intense erosion rates (∼ 4,020 t km⁻² yr⁻¹, ∼ 1.5 mm yr⁻¹ ). Similarly, sands show a slight increase of +4 [+1, +7] × 10⁶ t yr⁻¹ after accounting for Andean lateral inputs, resulting in a total flux increase of +14 [+3, +24] × 10⁶ t yr⁻¹. Floodplain drainage and sediment resuspension gain importance in this segment.

Contaya Arch – Requena (Box 3, Fig. 2b). The river shifts northeastward into the extensive Marañón Foredeep, a tectonic funnel terminating at the Iquitos Arch. In April, water storage reaches 9.4 [7.7, 11.1] km³ . The maximum flood duration reaches 125 d. The flood impacts both fine sediments and sands (Fig. 11), however with different processes. Fine sediment deposition (−16 [−23, −10] × 106  t yr⁻¹, 7 % of the load) is driven by prolonged floodplain residence time, while sand sedimentation (−11 [−13, −9] × 10⁶ t yr⁻¹) results from energy dissipation when the floodplain is active. This process, never documented before in Amazonian rivers, captures approximately 14 % [10 %, 20 %] of the sand flux.

Overall basin budget. Integrating these three compartments, the Ucayali floodplain traps 165 [138, 192] × 10⁶ t yr⁻¹ of the 455 [410, 500] × 10⁶ t yr⁻¹ supplied by the Andes (∼ 36 %), with sand accounting for 65 % of total deposition. The downstream export to the Amazon River reaches 290 [235, 345] × 10⁶ t yr⁻¹ (26 % sand), identifying the Ucayali as the dominant sediment source among Andean foreland tributaries of the Amazon, exceeding the Marañón (∼ 144–173 × 10⁶ t yr⁻¹; Armijos et al., 2013) and the Beni (∼ 133 × 10⁶ t yr⁻¹; Vauchel et al., 2017). It also exhibits some of the highest specific erosion and deposition rates in the foreland. Net fine sediment loss across the basin amounts to 21 % of the Andean input, driven primarily by the dampening of Andean flash floods in Compartment 1. Floodplain recycling during the recession limb contributes approximately 22 % of the total suspended load at the basin outlet, acting as a secondary sediment source that partially offsets net deposition and sustains sediment export during the low-water season. Peak floodplain water storage reaches 19.1 [15.3, 22.9] km³ in March, equivalent to ∼ 38 % of total basin discharge for that month.

6.4.2 Long-term trends

No significant trends were observed in flood duration, floodplain water volume, discharge extremes, or sediment fluxes in sub-basins 21, 22, and 23 for the simulation period. This suggests that the effects of global and regional changes at the Ucayali Basin scale are still negligible, contrasting with studies of the broader Amazon Basin (e.g. Marengo, 2004; Gloor et al., 2013; Li et al., 2020; Fleischmann et al., 2023). However, most of these studies focus on a very limited number of sites with long-term data. This highlights the relevance of the integrated approach proposed in this study for generating new long-term analyses of regional trends and providing a more detailed perspective on the Amazon's response to global changes. Espinoza et al. (2009) is the only study on Ucayali discharge evolution, based on inconsistent 1996–2005 data, now curated here. Figueiredo et al. (2024) and Polasky et al. (2025) identified mixed precipitation trends in the Ucayali Basin, respectively for the 1982–2023 and 1980–2022 period, which may explain the lack of clear trends here. Nevertheless, Guimberteau et al. (2013, 2017) project increased precipitation in the western Amazon by the century's end, potentially leading to greater inundations, sediment deposition, and floodplain recycling, which may alter river mobility, particularly during extreme events. This could impact bedload transport, leading to bed aggradation, which would facilitate overflow and then intensify exchanges.

6.4.3 Implications for sediment routing in large floodplain river systems

The floodplain-controlled sand sedimentation process identified here, whereby inundation reduces main-channel transport capacity and captures ∼ 14 % of the sand flux at peak discharge, generating a secondary sediment peak during recession, has implications beyond the Ucayali Basin. Large tropical river systems draining tectonically active forelands onto extensive alluvial plains are likely to exhibit similar morphological conditions, suggesting that floodplain-driven decoupling of sand flux from water discharge may occur more broadly. The integrated approach developed here provides a framework to identify and quantify this process in other poorly monitored basins.

7 Conclusion and perspectives

Three core findings emerge from this study.

First, floodplain inundation exerts a dual control on sediment dynamics: it captures ∼ 14 % [10 %, 20 %] of the sand flux at peak discharge by reducing main-channel transport capacity, while recycling processes during recession contribute 22 % of the total suspended load at the basin outlet. This floodplain-driven decoupling of sand flux from water discharge, previously undocumented in large Amazonian rivers, suggests that conventional approaches assuming a stationary relationship between water discharge and sand flux may overlook a first-order process in large floodplain river systems.

Second, the integrated framework developed here enables the first long-term hydro-sedimentary budgets for the Upper Amazon at sub-basin scale. The Andean Ucayali exports 455 [410, 500] × 10⁶ t yr⁻¹ of suspended sediment (40 % sand), of which 36 % is trapped within the floodplain, predominantly as sand (65 % of total deposition). Floodplain recycling contributes 22 % of the outlet load, and the Ucayali ultimately delivers 290 [235, 345] × 10⁶ t yr⁻¹ to the Amazon River (26 % sand), identifying the Ucayali as the dominant sediment source among Andean foreland tributaries of the Amazon. These estimates are stable over the 20–37-year record, with no significant long-term trends detected.

Third, in large floodplain rivers, the accuracy of sediment budgets is controlled by the sustained availability of multi-variable field observations at super stations, not primarily by model structure or parameter calibration. Despite substantial equifinality among secondary floodplain parameters, flux, trapping and recycling fractions remain stable across all behavioural GLUE simulations, demonstrating that long-term critical zone observatories such as HyBAm provide the irreducible constraint that no model refinement can substitute.

These findings highlight the need to rethink hydrological network operations through the concepts of super and virtual stations, combining long-term, high-quality in situ data with satellite measurements to constrain this integrated approach. Synchronizing targeted hydro-sedimentary gauging campaigns with satellite missions would significantly enhance model performance.

Moreover, the study underscores the importance of regionally focused and long-term analyses over broad, continental-scale assessments, which often introduce high uncertainties and may misguide sustainable development and mitigation strategies.

Looking ahead, this framework can be extended to other Amazonian basins facing similar hydro-sedimentary monitoring challenges, through targeted calibration campaigns and model refinements. The progressive deployment of satellite missions, particularly SWOT and Sentinel-3 for water levels and Sentinel-2 for high-resolution surface reflectance, will further enhance its spatial applicability by providing key constraints to calibrate floodplain activation thresholds and fine sediment concentrations in ungauged reaches.

The assimilation of satellite-derived fine sediment concentrations would significantly improve floodplain modelling by defining Floodplain Hydrological Response Units based on morphology (lakes, channels, etc.). By reconciling fine sediment mass balances (i.e., adjusting discrepancies between modelled and observed satellite data), this method would provide a proxy for water residence time in floodplains. This key variable is critical for constraining biogeochemical processes and could support the development of coupled bio-hydro-geochemical models, improving our understanding of the Amazon's role in global biogeochemical cycles and its response to environmental changes.

In conclusion, this study demonstrates that floodplains are not passive components of large river systems but active regulators of sediment transport, storage, and recycling. By integrating remote sensing, modelling, and long-term field observations, it provides a robust and scalable framework to quantify these processes, with implications for understanding biogeochemical cycles and managing large river basins under environmental change.

Appendix A: List of notations

Table A1 Main subscripts and superscripts.

Table A2 Main variables.