Flocculation Patterns Related to Intra-Annual Hydrodynamics Variability in the Lower Grijalva-Usumacinta System

Abstract

Particle aggregation modifies sediment dynamics, which is a determining factor for morphodynamic and ecological processes in deltaic plains. Here, we investigated the link between intra-annual hydrodynamics variability and flocculation in the Grijalva-Usumacinta system. Monthly (2016–2017) and seasonal (2021–2022) river data was processed using analytical methods and the simplified sonar equation. Flocs were reformed and characterized in the laboratory, validating the in situ settling velocities (0.5–3.8 mm/s) and the existence of large low-density macro-flocs (>300 μm). We verified that flocculation prevailed, exhibiting seasonal patterns; (1) the highest aggregation rates matched the increase in total suspended solids at rising-flow (>100 mg/L), (2) periods of high-flow showed stable aggregation rates, and (3) an influence of marine conditions occurred at low-flow. Particulate phosphorous and organic fraction showed seasonal patterns linked to flocculation. Due to damming, the shear rates varied slightly (7–11 L/s) in the Grijalva, leading to high flocculation intensities affecting the diffusivity ratio. In the Usumacinta, aggregation was limited by shear rates that normally exceed 15 1/s. We found seasonal Rouse parameters representative of sediment dynamics.

1. Introduction

Large fluvial systems drain the continental surface, gathering erosion products and anthropogenic agents, and form huge deltas that result from transport of sand, but mainly silt and clay [1,2,3]. These fine cohesive sediments (mud) host pollutants and bind important portions of nitrogen (N) and phosphorus (P) and organic carbon (OC) [4,5,6,7,8,9]. Because particles smaller than 32 μm tend to form aggregates (or flocs) [10,11,12], freshwater flocculation likely dictates sediment dynamics in large lowland rivers [13,14,15,16], and thereby, the dispersed particle approach does not appropriately represent the suspended load [17]. Likewise, tropical fluvial systems (30° N–30° S), show a marked intra-annual variability in flow (hydrological seasonality) that causes annual hysteresis patterns in the suspended load [18,19,20]. Therefore, when addressing morphodynamical processes and ecosystem health in floodplain and coastal areas, it is necessary to assess particle aggregation in river reaches and its link with hydrological seasonality [21,22,23]. This is not a straightforward task, since flocculation, even in freshwater environments, is a multifactorial process induced and affected by turbulence, organic matter content, sediments chemical composition, among others [24,25,26,27].

Ecological processes in rivers are linked to hydrodynamic forces at different time scales. For instance, intra-annual variability in flows (monthly scale) mainly influences lotic ecosystem structure and function at the community level [28,29]. On the other hand, hydropower is endorsed as a sustainable energy source; however, because of impoundments, the hydrodynamic regime can be drastically altered, affecting the sediment load and associated nutrient dynamics in the reservoir area and downstream to the plain [30,31,32]. Indeed, anthropogenic stressors such as changes in land use and damming may be even more critical for morphodynamic and ecological processes than seasonality [33]. The Grijalva-Usumacinta is a tropical fluvial system, regulated in the case of the Grijalva; therefore, this study addresses the implications of hydrological seasonality and damming on sediment dynamics in the plain.

We investigated the link between hydrodynamic forces and flocculation patterns with the intra-annual variability in flow at the lower reaches of the Grijalva and Usumacinta rivers, which combined represent the largest in Mesoamerica and the second input of freshwater into the Gulf of Mexico. We focused on flocculation implications for Rouse parameters and particulate phosphorous and organic matter. For this aim, shear velocity and dissipation-related parameters were derived from velocity data recorded monthly from May 2016 to June 2017. Flocs extracted from the river were reformed into a rotating annular flume, and their sizes and settling velocities were obtained with particle tracking velocimetry (PTV). In order to verify and calibrate the results derived from monthly data (2016–2017), vertical profiles of total suspended solids (TSS), organic matter content or organic fraction (OF), and particulate phosphorus (PP) were obtained seasonally (2021–2022). An empirical solution to the sonar equation was applied to derive the instantaneous TSS from backscatter signals, the mass conservation equation was then used to solve the in situ settling velocity. It is not the first time that similar methods have been applied, [34,35,36,37]; however, we used multiple approaches here, combining theoretical models, PTV, and the sonar equation to review the Grijalva-Usumacinta system. There are indeed few studies like this that complement field data with laboratory experiments to address flocculation in turbulent freshwater rivers [14,15,16,22].

Our working hypothesis assumed that the flocculation regime prevails, dictating the sediment transport in the lower Grijalva-Usumacinta. We addressed three main questions: (1) What are the patterns that characterize the flocculation regime during the different hydrological seasons? (2) How does flow regulation affect the suspended load dynamics in the lower Grijalva? (3) What implications does flocculation have on particulate nutrient dynamics? The main objectives were: (1) to verify that hydrodynamic forces follow annual hysteresis patterns similar to those of the suspended sediment flux (SSF); (2) to identify and explain the expected contrasts between the suspended load patterns in the Grijalva and Usumacinta rivers; and (3) to confirm that flocculation processes are mainly dictated by turbulence intensity and thus simplify the assessment of key parameters to model the suspended load dynamics in the plain.

2. Materials and Methods

2.1. Study Area

The Grijalva-Usumacinta is the largest tropical river-system in North America, representing 30% of México’s total runoff in just 4.7% of the country’s surface. Its watershed (121,025 km²) is located mostly in Guatemala and southeastern México and a portion in Belize (Usumacinta basin) [38]. This system drains the “Lacandona” forest (Chiapas) and part of the Petén (Guatemala). Headwaters are located in the tropical forest of Guatemala, and both rivers run through Chiapas and Tabasco in México to their confluence at “Tres Brazos” in the Tabasco plain before flowing into the Gulf of México (Figure 1). The Grijalva-Usumacinta delta is the 14th largest on Earth (21,000 km²) and part of the Gulf-Caribbean region of México that concentrates more than 65% of the forest reserves of the coastal plain, 75% of the coastal wetlands, and almost 50% of the total shrimp fisheries [39,40].

Mean annual temperature and average precipitation are 26° C and 3000 mm/year, respectively, but seasonality is notable, exhibiting a dry season (March–May), a rainy season (June–September), and a winter storms period (cold fronts) (October–February) locally called “Nortes”. The highest flows occur during the cold fronts period, reaching up to 8000 m³/s [40]. These high flows produce recurrent floods that cause economic damage and threaten the safety of the population [41,42]. For instance, just in November 2020 a serious flood occurred in the Tabasco plain, including Villahermosa, the most populous city in the Grijalva-Usumacinta basin (~1 million inhabitants).

2.2. Montlhy and Seasonal River Data

Sampling sites were located on the river main channel, at Los Ídolos and Chichicastle, Grijalva and Usumacinta, respectively (Figure 1a). Monthly samplings were carried out from May 2016 to June 2017, which comprised water velocities (ADCP data) and depth-integrated water samples (8 L) used to derive the depth-averaged TSS, OF and total phosphorus (TP) at 3 columns. In the case of Grijalva, samplings from July to February were carried out at S11 (4.5 km downstream from the Chilapa River entrance), whereas, due to the entry of the salt wedge, samplings from March to June (dry season) were conducted at S12 (upstream from the confluence with the Chilapa and the El Viento lagoon) (Figure 1b). Salt wedge did not penetrate as far as Chichicastle; thus, we sampled the Usumacinta at S21 every time (24 km upstream from Tres Brazos) (Figure 1c).

Within monthly sampling, seasonal data (assembled samples and TSS profiles) were obtained in the high-flow stage (08 November 2016: Grijalva, 01 December 2016: Usumacinta), and in the low-flow stage (dry season) (04 April 2017: both rivers. Assembled samples (50 L), they were collected by combining samples (2 L) from different depths at the channel main zone (0.2, 0.4, 0.6, and 0.8 h, where h = local depth) (Figure 1d–f). These were stored into plastic containers (10 L) and then transported to the laboratory (in coolers with ice). To complement the database, seasonal samplings were also launched on 19 November 2021 (high-flow), 04 April 2022 (low-flow), and 29 August 2022 (rising-flow), obtaining ADCP data and TSS, FO, and PP profiles (Figure 1d–f).

An Acoustic Doppler Current Profiler (ADCP: River Pro) tethered to a moving boat was used in bottom track mode to record flow velocities, conducting at least four transects for each measurement. Based on the flow conditions found by processing the monthly data, sampling columns for the 2021–2022 stretch were located in the most representative regions according to the hydrodynamic forces across the channel (Figure 1d–f). An isokinetic horizontal tube (Water trap), which is immediate-closing made of stainless steel (1.2 L, 2.5 kg of ballast) [43] (Lotsari et al., 2014), was used to obtain water samples (in triplicate) at five to six different depths for each column (Figure 1d–f). For each sample, TSS (concentration) and OF (organic matter fraction in suspended solids) were obtained applying the gravimetric method [44], using fiberglass filters (diameter = 4.7 cm) with a particles retention limit of 0.70 μm. Samples from different depths (lower, middle, and upper) at three columns were selected to determine (TP) and total dissolved phosphorus (TDP) (by persulfate digestion method), and soluble reactive phosphorus (SRP) (by ascorbic acid method) [44]. Particulate phosphorous (PP) was then calculated as PP = TP − TDP, [6].

2.3. Hydrological Regimen

Watersheds for the Grijalva and Usumacinta are similar in terms of geomorphology and hydrography [45]; however, these rivers run under different land uses and show contrasting hydrodynamic regimes. Since the mid-20th century, land use change has followed different routes in the middle and lower basins of these rivers [46,47]. Recent reports indicate that these contrasts are relevant to nutrient dynamics, affecting ecological processes in the plain [48,49,50,51]. Furthermore, the Grijalva (660 km in length) is regulated by four impoundments (hydroelectric dams) on its main channel (La Angostura, Chicoasén, Malpaso and Peñitas) [52,53]. The Usumacinta (1200 km in length), in contrast, is unregulated except for a dam on the tributary Chixoy River (Guatemala) in the upper basin.

In the case of Grijalva, seasonality was technically annulled by 1974 in the reaches upstream from Villahermosa (Carrizal River), that is, after La Angostura dam (the second) began to operate [52,53] (Figure 2a). At Los Ídolos, the Grijalva (main channel) has received the contributions from the unregulated rivers Pichucalco, La Sierra, and Chilapa (Tulijá and Macuspana rivers) (Figure 1). Therefore, flows measured at Los Ídolos (sampling site) are expected to reflect hydrological seasonality (Figure 2b,c). Boca del Cerro is located at Tenosique, Tabasco (80 km upstream from S21), being the nearest gauging station from Chichicastle (sampling site) (Figure 1). Because the Usumacinta is unregulated, gauged flows at Boca del Cerro adequately represent the hydrological regime in the main channel showing seasonality (Figure 2d).

2.4. Floc Characterization

To reform flocs in the laboratory, a rotating annular flume (1.30 m in diameter) was filled up to 5 cm in height with the water from the assembled sample, and a high shear velocity (>15 cm/s) was then generated during 60 min in order to induce disruption. After disruption, flocs were reformed under a shear velocity similar to that of field conditions (Grijalva: 4.5 cm/s, Usumacinta: 7.5 cm/s), which was maintained during 300 min. Floc sizes were measured every 15 min by capturing 90 double-pulsed images at each time interval. We assumed that floc size distribution reached at the end of the experiment (t = 300 min) corresponds to the floc-equilibrium state [54,55]. Once size distribution was obtained, the PTV array was moved to the settling column (1.20 m in height), and the annular flume was then restarted, generating again the reforming shear velocity for another 120 min. After 120 min elapsed, a ~60 mL water sample was taken from the flume and poured into the settling column, where floc sizes and their settling velocities were measured. Sets of 50 double-pulsed images were captured at the moment the sediment cloud was observed to reach the illumination sheet located at 80 cm deep. PTV algorithms applied to filter images (based on grey scale), and determine floc sizes and displacements, were derived from those used in similar studies, [56,57].

Floc size (equivalent sphere) was calculated as D_f = (4A_d/π)^1/2, where A_d = floc area according to image scale. Because they were derived from calibration, 1 pixel represented 47 μm; in the data analysis, we neglected flocs with D_f < 75 μm, that is, when A_d < 3 pixel. Floc settling velocity was calculated as ω_f = Δ_z/t_s, where Δ_z = displacement in the sinking direction (in mm), and t_s = distance between pulses in seconds (0.12–0.14). Displacement Δ_z was calculated from consecutive images, neglecting particles in resuspension (Δ_z < 0), and with a marked lateral displacement Δ_x > Δ_z, (Δ_x = displacement across the settling column). Further, flocs with a size variation lower than 30% between consecutive images, and ω_f in the range of 0.035 to 6.5 mm/s were considered reliable results.

2.5. River Hydrodynamics

Because turbulence is consistently three-dimensional (3D) in large rivers (due to transverse currents) [58], to calculate τ_r, we used τ_r = ρ(u′w′² + v′w′²)^1/2, where u′, v′, w′ denote the fluctuations in the velocity vector components, and u, v, w were streamwise, transverse, and vertical directions, respectively [59]. Shear velocity was then solved as u_∗ = (τ_r/ρ)^1/2 (ρ = water density). Likewise, u_* is commonly assessed by fitting the velocity profile to the Logarithmic Law, [36,60,61,62], estimating simultaneously the equivalent roughness, k_s = 30 y₀ (where y₀ = hydraulic roughness). Here u_* values derived from 3D-turbulent stress were compared/calibrated against the Logarithmic Law and then used to calculate the energy dissipation rate, ε = u_∗³/κy, where y = height above the bottom and κ is the von Karman constant (0.41) [63]. In this case, we calculated ε at y = 0.368h (the theoretical height where the mean streamwise velocity occurs). When quantifying the control exerted by turbulence on aggregation and/or disruption rates, two parameters derived from ε are used: the average shear rate, G, and the Kolmogorov microscale, η [54,64,65,66,67]:

$$G={\left(\frac{\epsilon }{\nu }\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}; \eta ={\left(\frac{{\nu }^3}{\epsilon }\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$4$}\right.}$$

(1)

where ν = cinematic viscosity. On the other hand, under varying turbulence conditions (changing hydrodynamic forces), τ_r is better represented either (1) by 3D-Turbulent kinetic energy (TKE): τ_r = 0.19 × ρ(u′² +v′² + w′²)/2 [59,61] or (2) by a modification to TKE: τ_r = 0.9 × ρw′² [68]. These expressions were used to calculate (calibrate) τ_r for the dry season (low-flow).

2.6. TSS Derived from Fluid Corrected Backscatter FCB

ADCPs can be used to directly measure SSF, thus improving the monitoring of river systems [69,70]. Furthermore, as in this study, knowing the instantaneous TSS allows the calculation of the Reynolds flux (sediment flux due to turbulence) [34]. We applied an empirical approach to derive TSS from FCB, following the simplifications considered by Venditti et al. (2016) in the Fraser River, United States, where concentrations and flows were similar to those measured in the Grijalva and Usumacinta rivers. Similar simplified procedures have been carried out on other large lowland rivers, [71,72,73]. Under this approach, the measured echo intensity, EI (in count), is converted to a backscatter signal, RL (in decibel), RL = s_f×EI (s_f = 0.43 is the scaling factor of the ADCP). The sonar equation is then applied to calculate FCB, that is, to correct RL for two-way transmission losses associated with fluid attenuation: FCB = RL + 20 log₁₀(r) + 2α_f, where r is the distance along the acoustic beam. The water absorption coefficient (under non-saline conditions) α_f = 8.686 × 3.38 × 10⁻⁶ (f²/f_T), where f_T = 21.9 × 10^{6−1520/(273−T)} (f = 1228 Hz is the frequency of the ADCP, and T is the water temperature in degrees Celsius). We assumed the two-way transmission losses associated with the suspended load were negligible because ambient TSS were low (<500 mg/L) [37], thus TSS can be derived from:

$$TSS={10}^{\left(A\times FCB+b\right)}$$

(2)

where A and b, are the coefficients (slope and interception) of the least squares regression. It is worth mentioning3 that this simplified application of the sonar equation ignores the effect of particle size distribution on backscattering and attenuation, [74,75].

2.7. Effective Settling Velocity

Under equilibrium conditions between upward turbulent diffusion and downward settling, the one-dimensional form of the mass conservation equation holds [76,77], hence the sediment flux due to turbulence (Reynolds flux) is equal to the effective settling velocity multiplied by the mean concentration:

$$\omega S=\overline{w\prime s\prime }$$

(3)

where ω = effective settling velocity; S = time-averaged concentration (mean TSS); and w′s′ = Reynolds flux (w′, s′ are the fluctuations in vertical velocity and concentration, respectively). Since the fluctuations, u′, v′, and w′ are known and the instantaneous concentrations are calculated by Equation (2), ω is then solved through Equation (3) [16]. Further, coefficients of momentum and sediment transport due to turbulence, ν_t and ν_s, respectively, can be calculated (ν_t = −u′w′/(du/dy); ν_s = s′w′/(ds/dy)), obtaining the diffusivity ratio (inverse of Schmidt Number), β = ν_t/ν_s [34,78].

2.8. Flocculation Intensity

Data obtained in the settling column (bin-averaged D_f and ω_f) was fitted to ω_f = AD_f^N, assessing the fractal dimension as F = 1 + N, and the representative floc size, D_r, by the effective settling velocity (ω), [35,54,79,80,81]. Khelifa and Hill (2006) introduced a model for floc-settling velocity (hereafter KH06) where the variation of density with D_f is addressed, and the primary particles size, d, is a predominant variable. Thus, to infer d, we also fitted the bin-averaged D_f and ω_f to KH06.

The suspension number (Rouse number), Z = ω/βκu_∗, indicates the concentration gradient, thereby the flocculation intensity [17]. For instance, in large lowland rivers, TSS gradients characterized by Z > 0.10 are presumably associated with high aggregation rates, where fine particles are transported as macro-flocs with ω > 2.0 mm/s [13,36]. In this case, for 3 (Grijalva) and 4 (Usumacinta) columns across the channel (Figure 1), we partitioned the analysis region in six bands of equal height (~0.125 h), obtaining (through Equation (2)) the mean S for each band at each column. Local Z were then solved by fitting S to the Rouse-Vanoni model [82,83]

$$S=S_a{\left(\frac{h-y}{y}\times \frac{a}{h-a}\right)}^Z$$

(4)

where S_a = reference concentration at y = a (in this case, a = height at the lower band ~0.15 h). Likewise, it is normally assumed that aggregation rate is limited (controlled) by turbulence intensity (in terms of G), leading to stable flocculation state where D_r/η < 1.0, [54,65,66]. Therefore, in contrast, a ratio D_r/η > 1.0 can be considered associated with a high aggregation rate characterized by low-density macro-flocs. Here, we take ω > 2.0 mm/s as a practical criterion indicating high flocculation intensity, verifying that it indeed corresponds to Z > 0.10, simultaneously assessing whether D_r/η > 1.0 is a common condition in the Grijalva and Usumacinta rivers.

3. Results

3.1. Intra-Annual Variability in Suspended Load

Flows measured from May 2016 to June 2017 were consistent with the seasonality observed at the gauging stations (Figure 1), with the lowest and highest flows matching the dry season and winter storms period, respectively (Figure 3). If it is assumed that measured data, q_i, reasonably represents the monthly average flow, then stages of rising, high, falling, and low flow can be inferred (Figure 3). If falling-flow stage comprises the winter storms period, then the flow measured on February, 2017 (Grijalva) was likely related to a cold front event. The mean annual flow (from May 2016 to June 2017), Q, was 634 m³/s and 1175 m³/s for the Grijalva and Usumacinta, respectively, with the Usumacinta contributing greater flows than the Grijalva, even for February (Figure 3a). Note that, for the Grijalva at low-flow (March, April, and May), q_i only takes into account the contribution of the main channel (Figure 1). Intra-annual variability in flow was evident, exhibiting relative flows, q_i/Q, ranging from 2.27 (October) to 0.06 (March) and from 1.75 (October) to 0.13 (April) for the Grijalva and Usumacinta, respectively (Figure 3b). As expected, flows (q_i) measured on November 2021 (high flow), April 2022 (low flow), and August 2022 (rising flow) followed the seasonality (Figure 3).

SSF showed annual hysteresis patterns in both rivers, although some relevant contrasts are worth addressing (SSF = q_i∙TSS). In the Grijalva, clockwise hysteresis in SSF occurred from the raising-flow to the falling-flow stage (from July to February), being interrupted at low-flow stage (q_i < 300 m³/s) (Figure 4a). In the Usumacinta, except for April, annual clockwise-hysteresis in SSF hold on (Figure 4b). Indeed, for the low-flow stage, SSF did not seem to match the hysteresis, but showed a wide range of variation (Figure 4a,b), likely derived from mean TSS that ranged from 18 mg/L (March 2017) to 70 mg/L (April 2022) and from 12 mg/L (April 2017) to 59 mg/L (April 2022) in the Grijalva and Usumacinta, respectively. Unlike SSF, OF did not show consistent patterns of intra-annual variation, except perhaps from the rising-flow to mid-falling stage (from June to December) (Figure 4c,d), where it can be observed that the higher TSS, the lower OF. However, even within this period, a marked decrease in OF occurred on November at both rivers (Figure 4c,d). From the end of the falling stage to low-flow stage (from January to May), OF varied randomly between 8 and 12% and 8 and 27% in the Grijalva and Usumacinta, respectively (Figure 4c,d). SSF, TSS, and OF measured in 2021 (November) and 2022 (April and August) were consistent in magnitude and with the trends of the 2016–2017 data (Figure 4).

3.2. TSS, OF, and PP Profiles

An increase in TSS with depth, z, was highlighted at both rivers for the rising-flow and high-flow, showing quietly similar gradients in relative TSS even though TSS measured on 2021 and 2022 were higher than those from 2016 and 2017 (Figure 5). In contrast, for the low-flow, except once at the Grijalva (April 2017), TSS ranged around the mean and did not show a consistent gradient (Figure 5a,b). It is noteworthy that the Grijalva showed higher increases in TSS with z (during high and rising-flow periods) than the Usumacinta (Figure 5), even though at these stages the mean TSS in the Usumacinta was 3.3 and 1.4 higher than in the Grijalva. In the Grijalva, the mean TSS was barely around 30 mg/L, but the TSS at z ~ 0.850 h (lower region) reached up 2.5 and 3.3 times greater TSS than at z ~ 0.20 h (near the surface region) for the high-flow periods in 2016 and 2021, respectively. This ratio even increased to 7.5 for the rising-flow when the mean TSS reached up 105 mg/L. In the Usumacinta, on the other hand, TSS at z ~ 0.85 h was 1.4 and 1.5 times higher than at z ~ 0.20h, for the high-flow 2016 and 2021, respectively; ratios similar to that obtained for the raising-flow (1.30 for mean TSS = 157 mg/L) (Figure 5).

In the Grijalva, for the high-flow and rising-flow, OF decreased with z, falling from 0.20 to 0.14, and from 0.16 to 0.11, respectively. For the low-flow, OF appears to be uniform in the vertical (0.10 in average) (Figure 6). In the Usumacinta, OF profiles did not present different patterns, and the averages were even similar (0.10, 0.09, and 0.09 for high, low, and rising flow, respectively) (Figure 6). In general, PP profiles were consistent with TSS at both rivers, increasing with z at the high-flow and rising-flow and ranging around the average for the low-flow (Figure 6). However, the change in PP with z was lower than that of TSS. In the Grijalva, PP at the lower region (z/h ~ 0.80) was 1.2 and 1.40 times the mean PP for the high-flow and rising-flow, respectively (Figure 6). In the Usumacinta, for the high-flow, the PP at the lower region was 1.10 times the mean PP, and it even decreased to 0.90 for the rising-flow (Figure 6). These patterns were related to the amount of P in the suspended solids, PP/TSS, which decreased with z at both rivers. In the Grijalva, PP/TSS of 1.5 and 1.7 mg/g in the upper region decreased to 0.8 and 0.6 mg/g in the lower region for the high-flow and rising flow, respectively (Figure 6). In the Usumacinta, PP/TSS decreased slightly with z, going from 0.7 mg/g in the upper region to 0.6 and 0.4 mg/g in the lower region for high-flow and rising-flow, respectively (Figure 6). Except for the low-flow, the mean PP/TSS were higher in the Grijalva than in the Usumacinta (Grijalva: 1.06 (high-flow) and 1.16 (rising-flow) mg/g; Usumacinta: 0.69 (high flow) and 0.62 (rising-flow) mg/g. For the low-flow (at both rivers), PP/TSS ranged around the average in the vertical, with means 0.60 and 0.68 mg/g for Grijalva and Usumacinta, respectively (Figure 6).

Particulate phosphorus represented between 0.50 to 0.95 of TP and varied between rivers and flow stage. In the Grijalva, PP comprised (in average) 0.60 and 0.80 of TP for the high-flow and rising-flow, respectively. In the Usumacinta, mean PP/TP was 0.92 and 0.72 for the high-flow and rising-flow, respectively. The lowest PP/TP values were obtained during low-flow, with 0.53 and 0.56 for Grijalva (TSS = 70 mg/L) and Usumacinta (TSS = 59 mg/L), respectively. In the Grijalva, for the high-flow (TSS = 35 mg/L) and rising-flow (TSS = 109 mg/L), the PP/PT ratio increased with the mean TSS. In the Usumacinta, in contrast, the PP/PT ratio for the high-flow (TSS = 110 mg/L) was lower than that of rising-flow (TSS = 157 mg/L). The ratio between SRP and TP also varied with the flow stage, presenting mean SRP/TP of 3 and 4% (Grijalva) and 1 and 4 % (Usumacinta) for the high-flow and rising-flow, respectively. The highest ratios SRP/TP were obtained for the low-flow, with 6% at both rivers. In average, the SRP/TDP ratios were 7, 13, and 23 % (Grijalva) and 8, 16, and 17% (Usumacinta) for the high-flow, low-flow, and rising-flow, respectively.

3.3. Hydrodynamic Forces

From the rising to the falling-flow stage, the velocity profile fulfilled the Logarithmic Law (0.94 < R² < 0.99), with logarithmic regions increasing with q_i at both rivers (Figure 7a,b). Indeed, within this stretch, the Logarithmic law predicted technically the same u_* as the 3D turbulent stress (Figure 7c); however, there is some relevant contrast worth addressing. In the Grijalva, a velocity dip pattern highlighted and holds up even for the high-flow stage (Figure 7a), whereas, in the Usumacinta, the logarithmic region covered almost the entire vertical (Figure 7b). At low-flow, except for the gravity current conditions, the streamwise velocities poorly fits the Logarithmic Law in the Grijalva (R² < 0.50) (Figure 7a) but holds on in the Usumacinta (R² ~ 0.90) with logarithmic regions covering at least 80% of the column (Figure 7b). Gravity current events were observed on March (39 m³/s, q_i/Q = 0.06) and April 2017 (155 m³/s, q_i/Q = 0.13) for Grijalva (overflow) and Usumacinta (underflow), respectively (Figure 7a,b). Under these conditions, a reasonable estimation of u_* was obtained with the modified TKE for the overflow event (Grijalva–March 2017), and in contrast, the 3D-TKE gave better results for the underflow event (Usumacinta–April 2017) (Figure 7c).

The semi-empirical solution for the sonar equation worked reasonably well, yielding acceptable TSS predictions (Figure 7d). However, different regression coefficients, A and b, applied for each river, which also varied with TSS and the flow stage; that is, it was not found a suitable general solution but particular ones depending on each river conditions. For the Grijalva, A, b, and R² were 0.036, −0.987, and 0.77 (high-flow: Nov/2016 and Nov/2021); 0.042, −1.217, and 0.82 (low-flow: April 2017); 0.024, 0.421, and 0.83 (low-flow: April 2022); and 0.077, −3.181, and 0.84 (rising-flow: August 2022). For the Usumacinta, A, b, and R² were 0.013, 0.860, and 0.78 (high-flow: Nov/2016); 0.066, −2.839, and 0.85 (low flow: April 2017); 0.022, 0.490, and 0.85 (high-flow: Nov/2021); 0.014, 0.932, and 0.77 (low flow: April 2022); and 0.027, 0.279, and 0.79 (rising-flow: August 2022). Further, for those months in which field TSS profiles were unavailable, the set A and b (by river) that best predicted the mean TSS in each case was used. In the case of July/2016 (both rivers), solving the sonar equation combining TSS from Nov/2016, Nov/2021, and April 2022 gave better results, yielding the coefficients A, b, and R² as follows: 0.035, −0.892, and 0.74 (Grijalva) and 0.029, 0.102, and 0.90 (Usumacinta).

In general, the Usumacinta runs under a higher turbulence intensity than the Grijalva, with u_* reaching up to 0.095 m/s (rising-flow and high-flow stages), whereas, in the Grijalva, u_* did not exceed 0.065 m/s (high-flow) (Figure 7c). Turbulence magnitude (in terms of G) follows hysteresis patterns similar to those showed by SSF. That is, in the Grijalva, the clockwise hysteresis occurred from the rising-flow to the falling-flow (from July to February) (Figure 8a). In the Usumacinta, except for April 2017, G followed a consistent annual clockwise hysteresis cycle (Figure 8b). For the Grijalva at low-flow, G varied strongly reaching magnitudes even higher than at high-flow (Figure 8a). This is due to values of u_* not less than 0.037 m/s reaching up to 0.049 m/s under gravitational current events in combination with a channel depth at location S12 that is half that at location S11 (Figure 2). In the Grijalva, the hysteresis in hydrodynamic force was slightly with G ranging between 11 and 7 1/s, for August (rising-flow) and February (falling-flow), respectively (Figure 8a). In the Usumacinta, hysteresis was remarkable, leading to G ranging between 6 and 24 1/s, for May (low-flow) and September (rising-flow), respectively (Figure 8b). Due to these patterns, η ranged from 276 to 383 μm (mean: 318 μm) and from 188 to 368 μm (mean: 248 μm) for the Grijalva and Usumacinta, respectively (Figure 8c).

3.4. Flocculation Prevalence

TSS profiles predicted from FCB matched those in the field, except for April 2022 (both rivers), fulfilling the Rouse model (R²: 0.60–0.97) with Z indicating a flocculation with widely range of intensities (Z: 0.06–0.70). For April 2022 (both rivers), contrary to the field data, predicted TSS followed the Rouse model. For the overflow event (Usumacinta: April 2017), TSS profiles did not fit the Rouse model (R² < 0.30). For the underflow event (Grijalva: March 2017), field TSS profiles were unknown, but as in the case of April 2022, predicted TSS followed the Rouse model. It can be said that the applicability of the Rouse model was better as the flow fitted the logarithmic law. We considered monthly averages for Z and ω, omitting those cases in which the Rouse model was obviously not applicable, confirming in summary, two patterns. (1) the flocculation regime prevailed throughout the year, and (2) the aggregation rate in the Grijalva was consistently higher than in the Usumacinta (Figure 9).

In the Grijalva, Z increased markedly during the rising-flow stage, from about 0.10 (May and June) to almost 0.60 (August). After this abrupt increase, Z decreased to 0.30 (September), then varied slightly, remaining close to the annual average (0.27) for the high-flow stage (Figure 9a). The highest ω also occurred at August, reaching 2.6 mm/s, although unlike Z, ω gradually decreased to 1.8 mm/s (December) (Figure 9b) with values all higher than the annual average (1.6 mm/s). It can be argued that the intra-annual trends of Z and ω were similar during the rising-flow (from June to August), with the highest Z and ω matching the highest mean TSS (> 100 mg/L) (Figure 9c). Afterwards a period of stability occurred for the high-flow stage (from September to December) (Figure 9a,b). These stretches presented high flocculation intensities characterized by Z > 0.25 and ω = 2.0 mm/s, even though the mean TSS decreased to ~30 mg/L for the high-flow stage [Figure 9 (c)]. Whereas, for the falling and low-flow stages, lower aggregation rates occurred (ω: 0.90–1.60 mm/s), and the trends of Z and ω do not necessarily match each other (Figure 9a,b)].

In the Usumacinta, the flocculation regime was relatively stable with Z and ω ranging between 0.07 (January) and 0.15 (August) and between 0.80 (April) and 1.50 (September) mm/s, respectively (Figure 9d,e), which are indeed comparable to the annual averages Z = 0.10 and ω = 1.2 mm/s. The highest aggregation rates occurred during transition from rising-flow to high-flow (between August and September), corresponding with the highest mean TSS (>150 mg/L) (Figure 9f). After the increase in TSS promoted flocculation, a period of stability in Z (~0.09) with randomly varying ω (0.80–1.4 mm/s) occurred from the high-flow stage to the end of the low-flow stage (from October to June) (Figure 9d,e). Within this stretch, the mean TSS ranged between 15 and 60 mg/L (Figure 9f). In summary, the highest flocculation intensities can be presumed to be associated with the rising-flow stage where TSS increases markedly in both rivers. In contrast, for the other flow stages, the flocculation intensity cannot be associated solely with the intra-annual variability in TSS.

On average, Z and ω for the Grijalva were 2.8 and 1.4 times higher than those of the Usumacinta, respectively, and these corresponded to values of G for the Usumacinta 1.8 times higher than those of the Grijalva. Furthermore, in the case of Grijalva, low Z (~0.10) occurring from May to July were related to G ~ 10 L/s (similar to those of the high-flow stage), even though TSS varied widely in this stretch (20–100 mg/L) (Figure 9a). For the Usumacinta, it can be argued that Z ~ 0.10 occurring between falling-flow and low-flow stages (from January to May) are perhaps explained by the decrease in G below 12 L/s (Figure 9d). Note that, according to annual-averaged conditions, Z ~ 0.10 represents a low aggregation rate for the Grijalva but an average for the Usumacinta. Likewise, G ~ 10 L/s perhaps means high turbulence intensity for the Grijalva but low for the Usumacinta. On the other hand, contrary to what might be expected, in both rivers from the high-flow to the falling-flow stage, the trend of Z seemed to follow that of G; that is, Z was positively related to G rather than negatively.

Regarding the intra-annual variability in ω, two patterns were notable in both rivers. First, the higher aggregation rates (in terms of ω) resulted from the increase in both G and TSS that characterized the rising-flow stage. Outside this stretch, the trend of ω did not seem to follow that of G (Figure 9). Second, between August and October, conditions of high ω (>2.1 mm/s: Grijalva, >1.3 mm/s: Usumacinta) corresponded to low η (Grijalva < 290 μm, Usumacinta < 210 μm). These high ω were likely related to large macro-flocs, and thus D_r/η > 1.0 perhaps occurred under such conditions. Note that, as in the case of Z, the high ω ranking refers to annual average aggregation rates, respectively. The trends of ω and OF did not appear to be related (Figure 9c,f). However, following Nghiem et al. (2022), we estimated the fraction of sediment surface covered by organic matter, θ, obtaining θ ranging from 0.07 to 0.35 (0.24 in average) and from 0.02 to 0.26 (0.14 in average), for the Grijalva and Usumacinta, respectively. These θ can be classified as high compared to 0.07 ± 0.04 (global average in large rivers) [67] and, therefore, a promoting factor for particle aggregation. Thus, the higher OF (or θ) in the Grijalva than in the Usumacinta supports the higher aggregation rate in the Grijalva.

3.5. Diffusivity Ratio and Floc Size

Factor β proved to be a relevant parameter to adequately apply the Rouse model in these rivers. In both cases, β was negatively correlated with u_*/ω, resulting in β = 4(u_*/ω)^−0.80 (R² = 0.62, Grijalva) and β = 16.26(u_*/ω)^−0.91 (R² = 0.79, Usumacinta) (Figure 10a). In the mean, β in the Grijalva (0.32 ±0.12) was lower than in the Usumacinta (0.54 ±0.33), indicating a contrast between their sediment loads highlighted by Z(Grijalva)/Z(Usumacinta) > ω(Grijalva)/ω(Usumacinta). Data was contrasted to the empirical expression introduced by de Leeuw et al. (2020) for sand mixtures (compact grains) and extended by Lamb et al. (2020) for rivers under a flocculation regime, β = 16.82(u_*/ω)^−0.54C_f^0.3 (C_f = u_*²/U², U = mean streamwise velocity) (Figure 10b). Note that in this case, we did not calculate β by fractions, e.g., [17,67], but assumed that the suspended load is mainly composed of clays and fine silts (<16 μm), and therefore, flocs represent the mixture. Likewise, instead of the skin-friction velocity we used the mean shear velocity (u_*). For the Usumacinta, the data were acceptably predicted and well correlated (R² = 0.81), resulting in β(predicted)/β(Usumacinta) = 1.3 +/−0.2. For the Grijalva, on the contrary, the factor β was notoriously overestimated and poorly correlated (R² < 0.10), resulting in β(predicted)/β(Grijalva) = 2.8 +/−1.1.

Bin-averaged ω_f and D_f followed the empirical relationship ω_f α D_f^(F−1) (R² > 0.80) (Figure 10c,d), resulting in F reflecting the turbulence intensity in both rivers. Based on these fits, for the Grijalva, the mean F equaled 1.9 and 1.3, correspond to the flocs formed at the high-flow stage (November/2016) and the low-flow stage (April/2017), respectively. While for the Usumacinta, F equaled 2.4 and 2.2, which corresponded to the high-flow stage (December/216) and the low-flow stage (April/2017), respectively. Because ω and G were quite similar for months comprising the high-flow and low-flow stages, respectively (Figure 9), we extrapolated the empirical fit corresponding to the high-flow to estimate D_r for October, November, and December and that of the low-flow stage for March, April, and May. It resulted in D_r ranging from 625 to 725 μm (high-flow) and from 70 to 205 μm (low-flow) for the Grijalva and from 330 to 415 μm (high-flow) and from 125 to 495 μm (low-flow) for the Usumacinta. These D_r were in the range of macro-flocs (>300 μm), except for the low-flow stage in the Grijalva, resulting in D_r/η ratios (in average) 2.1 (high-flow) and 0.4 (low-flow) for the Grijalva and 1.8 (high-flow) and 1.6 (low-flow) for the Usumacinta.

Based on KH06, for the Grijalva (high-flow), d = 1.40 μm worked for D_f in the range of 75 to 250 μm (ω_f: 0.30–1.05 mm/s) while d = 8 μm was better suited for large macro-flocs (D_f: 500–700 μm, ω_f: 1.20–3.40 mm/s) (Figure 10c). In the case of the Usumacinta, for the high-flow and low-flow stages, the same d = 0.70 μm worked for D_f in the range of 75 to 270 μm (ω_f: 0.15–0.50 mm/s) (Figure 10d). For large macro-flocs, in contrast, the best fits resulted in d = 8 μm for the high-flow (D_f: 600–900 μm, ω_f: 3–4 mm/s) and d = 2.5 μm for the low-flow (D_f: 375–700 μm, ω_f: 0.97–2.55 mm/s) (Figure 10d). For the Grijalva (low-flow), a single d = 4.5 μm worked for almost the entire population (D_f: 85–580 μm, ω_f: 1.15–2.25 mm/s) (Figure 10c). In all cases, the best fits were obtained for a primary particles density, ρ_s = 2300 kg/m³, resulting in effective densities, ρ_e, ranging from 273 to 13 kg/m³ (high-flow) and from 264 to 20 kg/m³ (low-flow) for the Grijalva and ranging from 302 to 11 kg/m³ (high-flow) and from 209 to 10 kg/m³ (low-flow) for the Usumacinta (ρ_e = ρ_s − ρ_f, ρ_f = floc density).

4. Discussion

4.1. Study Findings

Here, it was verified that hydrodynamic forces do indeed follow annual hysteresis patterns, matching those of SSF. Remarkably different effective settling velocities characterize the sediment dynamics in these rivers, mainly explained by the intra-annual variation in G, also indicating different suspended load compositions. For instance, diffusivity ratios (β) were consistently lower than 1.0 and were negatively correlated to u_*/ω, but with different regression coefficients for each river (Figure 10a). These results confirmed that fine particles (in floc forms) indeed affect coefficients ν_t and ν_s [84], and furthermore, β depends on floc properties (D_f, ρ_e) and perhaps on suspended load composition. Regarding D_r, if indeed D_r is positively related to η, flocs formed in the Grijalva are expected to be larger than those in the Usumacinta, expecting a higher intra-annual variability in D_r for the Usumacinta than for the Grijalva. Further, enough evidence has been found here to support D_r/η > 1 as a common condition, except in the Grijalva for the low-flow stage. Therefore, it is recommended to characterize flocs from the rising and falling stages. In the absence of data, it seems reasonable to assume a D_r/η similar to that of the high-flow and low-flow stages for the rising and falling stages, respectively.

4.2. Seasonality and Flocculation Patterns

Flocculation patterns and their intra-annual trends were identified, resulting in the assessment of Rouse parameters useful for predicting SSF and even sediment dynamics when field TSS profiles are not available. From a practical perspective, three stretches highlighted for the Grijalva: (1) low-aggregation rate (February to June: falling and low-flow stages), (2) maximum flocculation intensity (August: at the peak of the rising-flow stage), and (3) stable-high flocculation (from September to December: end of rising-flow stage and the entire high-flow stage), which can be characterized by a mean ω equal to 1.2, 2.5, and 2.0 mm/s, respectively. For the Usumacinta, even though it presents a marked annual hysteresis in G, two aggregation periods seem sufficient: (1) controlled aggregation rate (from October to July: from the high-flow stage to the beginning of the rising stage) and (2) promoted aggregation (August and September: rising-flow stage), which can be represented by a mean ω equal to 1.1 and 1.5 mm/s, respectively. Even for rough SFF predictions, the mean annual ω (1.2 mm/s) is a reasonably proxy for the suspended sediment mechanics.

On the intra-annual variability in flocculation, there were two patterns worth addressing: (1) the positive relationship between Z and G occurring from the rising-flow stage to the falling-flow stage (both rivers) and (2) the high aggregation rates (Z > 0.20, ω > 2 mm/s) under low TSS (<40 mg/L) conditions in the Grijalva. From these, it follows that (1) the higher G, the lower Z does not necessarily hold for the intra-annual variation in the suspend load dynamics and (2) TSS is perhaps not the main factor promoting aggregation in turbulent freshwater environments. Likewise, when comparing the Grijalva and the Usumacinta, it reveals that (1) for G between 7 and 11 L/s, turbulence promotes aggregation leading to large macro-flocs and (2) a G > 15 L/s likely limits aggregation (Z < 0.15, ω < 1.5 mm/s), even for TSS > 150 mg/L. In other words, in the Grijalva–Usumacinta system, flocculation was mainly dictated by turbulence, leading to different patterns derived from two basic hydrodynamic conditions: (1) an optimum G range for aggregation (under which the Grijalva regularly runs) and (2) a limit G that controls aggregation (normally exceeded in the Usumacinta). This statement is consistent with the results reported by Li et al. (2020) for flocs from the Yangtze River (Three Gorges Reservoir), suggesting G > 16.5 as the limit where turbulence goes from promoting aggregation to regulating it.

Discriminating between controlled and high intensity flocculation states can be quite practical. For instance, when insufficient data is lacking, a ω ~ 1 mm/s (controlled flocculation) is commonly adopted when addressing large-scale morphodynamical processes involving sediment cohesiveness [3]. This is certainly realistic, but the predictions could be compromised under high aggregation rates. In the case of Grijalva, flocculation patterns such as D_f > 500 μm, ω > 2.0 mm/s, and F < 1.9 are unusual compared to average conditions in large rivers (D_f ~ 130, ω ~ 1.8, F > 2) [67], explaining the low correlation between the predicted β and the obtained β (Figure 10b). Large rivers around the world are vastly different in terms of flow rate, suspended load, chemical composition, organic matter content, and human impacts [85,86], but on the contrary, hydrodynamic conditions such as u_* and ε (G and η) are located in a very similar range of variation [17]. Thus, addressing the link between intra-annual variability in turbulence and flocculation processes is likely the most convenient option to improve the modelling of morphodynamic and ecological processes in floodplains and coastal areas. Here, we suggest verifying that there are indeed (1) a range of G promoting aggregation and (2) the limit G that prevents high aggregation rates.

4.3. Flow Regulation

It was argued that G ranging from 7 to 11 L/s caused high aggregation rates in the Grijalva. We hypothesize that the narrow range of variation in G in the Grijalva was likely related to the regulated flow in the middle reaches, causing the flattening observed in the G vs. q_i curve (Figure 8a). At the annual hydrographs (Gonzalez gauging station), it is evident how the hydrological regime has been altered, flattening the variability in flow rate (Figure 2a), with impacts reflected even downstream to Los Idolos. In other words, even with contributions of the unregulated rivers, La Sierra and Chilapa, that return seasonality to the lower Grijalva, the intra-annual variability in the hydrodynamic forces is reduced by the effect of impoundments, affecting sediment dynamics downstream to the floodplain. Damming rivers alters the hydrodynamic regime, further increasing the feasibility of particle aggregation in the reservoir area [87,88]. This not only affects the dam storage capacity, but also ecological processes [30,31,32]. Therefore, flocculation in the middle Grijalva is very likely, and its implications on sediment dynamics at the basin level must also be addressed.

4.4. Implications on Particulate Nutrients on SSF

We assume that the decrease in OF and PP/TSS with z resulted from the negative relationship between d and the amount of phosphorus and organic matter adsorbed by suspended sediments, as previously documented [5,8,26,89]. In other words, due to flocculation, D_f and d surely increased with z, causing OF and PP profiles that did not match that of TSS for the high-flow and rising-flow stages. This assumption is indeed consistent with the hypothetical variation in d (derived from KH06), suggesting that in this system, clays and silts (0.7–8 μm) formed large low-density macro-flocs (>300 μm). Likewise, also derived from KH06, the range of variation in d was presumably narrower for the low-flow than for the high-flow stage (Figure 10c,d). This result is consistent with the almost uniform TSS, OF, and PP profiles for the low-flow stage (both rivers). In light of these patterns, intra-annual variability in flocculation processes must not be ignored when addressing the particulate nutrient dynamics.

SSF and suspended load dynamics are commonly addressed under the wash load approach (neglecting flocculation) [18,19]. Assuming this mode of transport for the lower Grijalva-Usumacinta system leads to the assumption that fine sediments reach the mouth far out into the Gulf of Mexico. This hypothesis was strongly rejected by the evidence gathered in this study, highlighting relevant contrasts regarding the flocculation patterns in the Grijalva and Usumacinta rivers. Due to particle aggregation, the effective settling velocities can reach up two orders of magnitude to that of dispersed particles, reducing drastically the retention time and thus the transport distance. Therefore, under low-flow conditions, when the river transport capacity decreases allowing the salt wedge enters, fine sediment deposits can be produced along the floodplain and near the mouth. During low-flow periods, gravity currents may be common rather than the exception, leading to changing hydrodynamic conditions where flocculation prevails, but Logarithmic Law and Rouse’s model do not hold. That is, the almost uniform profiles of TSS, OF, PP/TSS obtained during low-flow are the result of gravity currents, and therefore, the analysis of sediment dynamics must take into account the periods of overflow and underflow [90].

Outside the low-flow stage, the suspended load presumably reaches the mouth, entering the Gulf of Mexico. However, in the light of this study, the mechanical properties of the suspended load are quite different before joining in Tres Brazos, characterized by a ω in the Grijalva that is twice that of the Usumacinta. Flocs resulting from this mixture are then transported towards the mouth, entering a decaying shear field [65]. Therefore, to improve the understanding of the ecological processes in the coastal zone, the flocculation processes at the confluence of the system must be addressed. On the other hand, as confirmed in this study, from rising-flow to falling-flow stages, the Rouse model is applicable to calculate the SSF.

4.5. Limitations, Strenghts, and Pending Tasks

A set of simplifications should be taken into account when analysing the findings of this study. (1) Primary particles size was not measured in the experiments, but rather inferred considering a floc settling velocity model [79]; (2) the resolution of the PTV array allowed only the measurement of flocs larger than 75 μm; and (3) no bottom samples were available. Therefore, the compositions of suspended and bed loads could not be compared. In contrast, the strengths of this study were that (1) the different flow stages that characterize the hydrological regime were addressed and therefore compared in terms of sediment dynamics, (2) the effective settling velocity was contrasted to settling velocities of real flocs (reformed in the annular flume), and (3) data processing was based on widely supported theoretical models, confirming their applicability in both rivers.

Not having found a unique solution for the sonar equation could complicate future estimates of TSS from FCB; that is, coefficients A and b may not be extrapolated to other locations even within the same river system. These results were expected due to intra-annual variability in river hydrodynamics and suspended load composition. For instance, in the Fraser River, United States [37] and in large rivers in Russia [63], A and b also varied with TSS and flow conditions. Therefore, testing a calibration process involving floc properties perhaps allows the generalization of the TSS prediction, for which one needs to know in addition to D_f and ω_f, d [74]. Nonetheless, we used the simplified approach as a practical option, deriving results consistent with both theoretical background and experiments.

Regarding the flocs characterization, there is no doubt that it remains to complement the findings by addressing ω_f vs. D_f for flocs from the rising-flow and falling-flow stages. For this, combining PTV with Digital Holography represents an interesting option, assessing D_f and the floc 3D-geometry and improving the PTV resolution [91]. Likewise, it is recommended to collect and analyze bottom samples and, for the suspended load, to address the elemental molar ratio Al/Si, the relative charge density, and the extracellular polymeric substances (EPS) [67]. These tasks should allow the correlation of the suspended load dynamics with the bottom transport, and then the derivation of a flocculation based-model, which would be useful to assess the morphodynamic and ecological processes in the Grijalva-Usumacinta plain.

5. Conclusions

The results of this study confirmed that the flocculation regime prevails in the Grijalva–Usumacinta system, with patterns linked mainly to intra-annual variability in hydrodynamic forces. This notably affects the dynamics of suspended load and associated nutrients, highlighting the following:

(1) The different aggregation rate (flocculation intensity) in the Grijalva and Usumacinta rivers is due to the range of variation in their G. In this case, G varying between 7 and 11 L/s and greater than 15 L/s presumably caused high aggregation rates in the Grijalva (Z > 0.25, ω > 2.0 mm/s) and controlled aggregation rates in the Usumacinta (Z ~ 0.10, ω ~ 1.20 mm/s), respectively.

(2) TSS, OF, and PP profiles are affected by the aggregation rate, presenting seasonal variation patterns and suggesting that, as in large lowland rivers, d increases with z, except during the low-flow stage.

(3) During the low-flow regime, the freshwater flow is blocked by the salt wedge, affecting the hydrodynamic conditions and thus flocculation patterns (TSS profiles), disabling even the applicability of the Rouse model. That is, estuary-like conditions occur even in the freshwater reaches.

(4) High aggregation rates, like those in the Grijalva, may not be considered when addressing SSF and morphodynamic processes based on common generalizations of flocculation processes. For instance, the factor β is evidently affected (decreased) when dealing with large low-density macro-flocs (ω > 2.0 mm/s, F < 2). For the Usumacinta, in practical terms, the controlled aggregation rate makes it reasonable to assume that ω ~ 1.2 mm/s remains relatively stable throughout the year.